NASA— CR- 192 349 


/y/ , ■ j - 7 3 0 
/// - 75 -C/Z 

/ *- 



MIXED CONVECTIVE/DYNAMIC ROLL VORTICES AND THEIR EFFECTS ON 
INITIAL WIND AND TEMPERATURE PROFILES 


Tracy Haack 
and 

Hampton N. Shirer 


Department of Meteorology 
The Pennsylvania State University 
University Park, PA 16802 


GO 



o 



(— ♦ 




i/> 

n 

<M 

n? 

lT 

1 


o 

m 

u 

in 

O' 

c 

.—4 


o 

o 


fSJ 

O 

m 

o 


^Present affiliation: 

Naval Oceanographic and Atmospheric 
Research Laboratory, 
Monterey, CA 93943-5006 


o 

iXJ 

U * 


QC < 


n *"• 

> k- * 

cn > 

a _j ^ ll -- 

UJ ^ -J c 

X o ^ 

M of 7 U- 

x n n v 

O Cl k 1 

t— « vl Cl 'V 

^ X k- V 

<> < lu t/» 

sf T' LL» 

m > ll 7> c? 

^ Q LL f- ■- 

O X. UJ -X c 

— LL aT »t 

| > c< LU > 

^ M k-i a *- 

o k LU y >* 

I O X Uj tn 

< LL: k I — C 

uo > C 

*T e O 73 3 

e C‘ 2* a. 

w o -X < w 


July 1991 



1 


ABSTRACT 

The onset and development of both dynamically and convectively forced boundary layer 
rolls are studied with linear and nonlinear analyses of a truncated spectral model of shallow 
Boussinesq flow. Emphasis is given here on the energetics of the dominant roll modes, on the 
magnitudes of the roll-induced modifications of the initial basic state wind and temperature 
profiles, and on the sensitivity of the linear stability results to the use of modified profiles as 
basic states. It is demonstrated that the roll circulations can produce substantial changes to 
the cross-roll component of the initial wind profile and that significant changes in orientation 
angle estimates can result from use of a roll-modified profile in the stability analysis. These 
results demonstrate that roll contributions must be removed from observed background wind 
profiles before using them to investigate the mechanisms underlying actual secondary flows 
in the boundary layer. 

The model is developed quite generally to accept arbitrary basic state wind profiles as 
dynamic forcing. An Ekman profile is chosen here merely to provide a means for easy 
comparison with other theoretical boundary layer studies; the ultimate application of the 
model is to study observed boundary layer profiles. Results of the analytic stability analysis 
are validated by comparing them with results from a larger linear model. For an appropriate 
Ekman depth, a complete set of transition curves is given in forcing parameter space for roll 
modes driven both thermally and dynamically. Preferred orientation angles, horizontal 
wavelengths and propagation frequencies, as well as energetics and wind profile 
modifications, are all shown to agree rather well with results from studies on Ekman layers as 
well as with studies on near-neutral and convective atmospheric boundary layers. 
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1. Introduction 

Observational evidence of stable secondary flows in the planetary boundary layer is 
abundant. Early scientists (Langmuir 1938) noted that long rows of seaweed floated parallel 
to the wind, and Woodcock (1942) observed that the soaring patterns of seagulls correlated 
with convective updrafts. More recently, satellite pictures have shown lines of sand dunes in 
the Sahara aligned with the prevailing winds and numerous examples of cloud streets during 
cold air outbreaks over water (Brown 1980). Atmospheric roll circulations consistently occur 
in boundary layers having a slightly unstable stratification and moderately strong winds. 
Typical wavelengths for these rolls range from two to eight km and are about three times the 
circulation depth; roll orientations are approximately 15* to the left of the geostrophic flow; 
phase speeds average one to two m/s; and propagation periods may be from 15 min to 2 hrs 
(Brown 1972; LeMone 1973). 

As reviewed extensively elsewhere (Brown 1980; Stensrud 1987a), roll circulations are 
driven by both convective and dynamic instability mechanisms. Convective instability 
produces Rayleigh/Benard circulations when a vertical potential temperature gradient exists 
at the surface. For wind speeds greater than a few meters per second, these thermal cells 
align in linear cloud bands that are nearly parallel to the direction of mean wind shear 
(Kuettner 1959, 1971). Studies of neutral atmospheres (Lilly 1966; Faller and Kaylor 1967) 
indicate that two dynamic mechanisms can also induce secondary flows. The inflection point 
instability mechanism generates roll circulations when sufficient energy is extracted from the 
sh ea r in the roll— perpendicular wind component. In contrast, rolls excited by the parallel 
instability mechanism require Coriolis turning to extract energy from the roll-parallel wind 
component; typically, the Coriolis conversion terms are much smaller in magnitude than are 
the other energy contributions, and so the parallel instability mechanism is believed to be of 
lesser importance in the atmosphere (Brown 1972; LeMone 1973, Brummer 1985; Chlond 
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1987). Brown (1980) notes that the convective and inflection point instability mechanisms 
appear to be sufficient to explain most geophysical lineal flow patterns. 

Each of these mechanisms has been studied extensively using a variety of theoretical 
approaches. Lilly (1966) and Brown (1970) investigated the dynamic mechanisms using 
linear models based on the partial differential equations for the neutral Ekman boundary 
layer. Etling (1971), Brown (1972), and Asai and Nakasuji (1973) are among those who 
utilized linear partial differential models to study the mixed convective/dynamic instabilities 
of an Ekman layer, while others, including Kuo (1963), Asai (1970a, 1970b, 1972), Kuettner 
(1971), Sun (1978), and Shirer (1980) considered the linear boundary layer responses to other 
sheared flows. In addition, high resolution nonlinear models have been developed for 
numerical study of both the Ekman layer (e.g. Faller and Kaylor 1966, 1967) and more 
general boundary layers (e.g. Faller and Kaylor 1969; Sommeria and LeMone 1978; Mason 
and Sykes 1980, 1982; Becker 1987; Chlond 1987; Etling and Raasch 1987). Results from the 
studies mentioned above have provided a basic understanding of: 

(i) the instability mechanisms responsible for roll development, 

(ii) the preferred roll characteristics, i.e. orientation angles and 
horizontal wavelengths associated with each mechanism, 

(iii) the resulting secondary flow patterns, and 

(iv) the profiles of the energetics terms, vertical transports and 
modifications to the basic state. 

A host of studies that incorporate measurements of observed boundary layer circulations (e.g. 
Kuettner 1971; LeMone 1973; LeMone and Pennell 1976; Sommeria and LeMone 1978; 
Weston 1980; Kelly 1984; Walter and Overland 1984; Becker 1987; Chlond 1987; Etling and 
Raasch 1987) have in general confirmed these theoretical results. 

Although the above studies have revealed a great deal about the fundamental properties 
of boundary layer rolls, they have provided only snapshots of the expected behavior in the 
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forcing and response parameter spaces. In particular, large, high-resolution models can not 
provide a very complete picture because they can be integrated only for a limited number of 
parameter values; moreover, in many cases the modeled rolls are not able to respond by 
changing their orientation angle because usually this angle must be specified in advance. In 
contrast, most linear models provide the response parameters as output, but they typically 
are developed for study of only one, usually the Ekman, profile and so their linear stability 
results can not be generalized easily to that expected for an atmospheric profile. Hence, a 
simple model capable of accepting arbitrary profiles is required for investigation of observed 
boundary layer flows. Existing attempts at comparing observations with model results tend 
to use the observed, roll-modified wind profiles as input rather than the (probably 
unobserved) pre-roll state whose instability actually led to the rolls themselves. The 
sensitivity of the linear analysis to use of such incorrect profiles must be investigated further; 
if significant sensitivity is found, then a simple means for estimating the probable roll 
modification is needed. Therefore, a modeling approach must be used that allows 
representation of the roll modes with a model large enough to capture the important dynamic 
and thermodynamic modes but still small enough to study analytically. 

Such a modeling approach is codified in the low-order spectral technique pioneered by 
Lorenz (1963) and discussed extensively in the book edited by Shirer (1987). In these 
nonlinear models, the dependent variables are represented by truncated Fourier expansions 
composed of temporally dependent amplitude coefficients and spatially dependent 
trigonometric basis functions. For boundary layer roll studies, the spatial characteristics are 
often given by a single horizontal harmonic and one (Shirer 1980, 1986) or two (Stensrud and 
Shirer 1988) vertical harmonics. The possibility of successfully using such severe truncations 
is suggested, for example, by ground observations (Kuettner 1959, 1971), examination of 
satellite images (Brown 1980), time-height cross sections of tower data (LeMone 1973), and 
cross sections given by aircraft data (LeMone and Pennell 1976; Briimmer 1985). Basing the 
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truncated models on the complete nonlinear Boussinesq equations allows study of both single 
and mixed instability mechanisms over a wide range of environmental wind shears and static 
stabilities (Shirer 1986). The onset of a roll mode is represented by a bifurcation from a 
motionless conductive state, and each critical forcing parameter value is given by a root of an 
analytically derived polynomial equation (Stensrud 1987a). Consideration of a large range of 
parameter space can therefore be performed quickly and efficiently. Moreover, because the 
basic state in these models is represented by a truncated Fourier series (Stensrud and Shirer 
1988), the spectral modeling approach is ideal for direct comparisons between model results 
and observations, as was done by Shirer and Briimmer (1986) and Stensrud and Shirer 
(1988). Although reasonable agreement was obtained between their results and observations 
taken during the 1981 KonTur experiment (Brummer 1985), the models were somewhat 
limited because the potential roll modification of the initial background wind could not be 
considered. 

In order to investigate the possible modification of the initial basic state by boundary 
layer rolls, we develop in sections 2, 3 and 4 a new nonlinear 14-coefficient spectral model of 
two-dimensional shallow Boussinesq flow that is forced both convectively and dynamically. 
In section 5 we qualitatively compare our results with those of previous theoretical and 
observational studies, and in section 6 we investigate the sensitivity of the stability results to 
the use of roll-modified profiles. In our study, we use an Ekman profile because it has been 
used in the vast majority of previous theoretical studies; however, the ultimate application of 
the model is to observed atmospheric profiles. 
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2. Model development 

To model boundary layer roll circulations, we use the shallow Boussinesq equations 
(Dutton and Fichtl 1969). Cloud streets typically form in the upward branches of these 
circulations, indicating the presence of two-dimensional roll patterns. However, we may 
neglect latent heating effects by assuming that the cloud area is small (Laufersweiler and 
Shirer 1989). The rolls are represented by perturbations superimposed on a time- 
independent, hydrostatic and horizontally moving basic state (Shirer 1980, 1986). Here, the 
initial cross— roll wind profile U(z) is approximated by a truncated Fourier series involving 
two vertical wavenumbers q and n: 

U(z) = Uo + Uisin(q irz/z,p) + U2 cos(q 7 rc/z r j-,) 

+ U 3 sin(n 7 rz/zrj.) + U 4 cos(n irz/zrj,) (2.1) 

where Zrj, specifies the height of the domain. The coefficients Ui are found from appropriate 
Fourier integrals, as described in section 5. As a first-order approximation, the initial 
temperature profile T 0 (z) varies linearly with height and is defined in two parts: 

T 0 (z) = Ti + Tf = (T sa - 7 ez) + (Tib - Tsa)(ziji— z)/z,p (2. 2) 

The temperature Ti represents the contribution owing to the environmental lapse rate 7 e and 
the temperature Tf is a vertically distributed surface forcing contribution based upon the 
difference between the lower boundary temperature Tib and the surface air temperature T sa . 

We assume that the perturbations possess a two-dimensional structure and so neglect 
all roll-parallel variations. Thus, the horizontal and vertical equations of motion may be 
combined into a single vorticity equation using the following form for the stream function: 
dip/dz = -u' and d^j dx. — w' . We also make the standard assumption (e.g. Brown 1970) 
that the Coriolis terms are small in magnitude and so do not contribute significantly to roll 
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development. Although Brown (1970) proposed a profile-modification mechanism that 
depended on the Coriolis parameter in his mean wind equations, we investigate here whether 
a more rapid adjustment mechanism is possible that depends only on the nonlinear coupling 
terms. Evidence for relatively rapid changes to observed cloud streets, implying rapid 
changes to the background profiles, is given by Brummer (1985). The above simplifications 
are strictly valid provided that the perturbations reach a steady state within a short time 
scale, two to four hours, and that only moderate supercritical forcing rates are considered in 
the temporal integrations. 

The horizontal domain is infinite and cyclically continuous at x = 0 and x = L, where L 
is the roll wavelength. Vertically the domain ranges from z = 0 to z = Zrp. For simplicity 
these boundaries are assumed to be rigid, stress-free and perfectly conducting. In 
dimensionless variables denoted by an asterisk, we have a domain defined by 0 < x* < 2x and 
0 < z* < 7T. 

Using the above assumptions, we derive a partial differential system containing 

J 

equations for the perturbation stream function ip and perturbation temperature T . We write 
the dimensionless equations representing boundary layer flow as (Stensrud 1987a; 
Haack-Hirschberg 1988) 

+ K' ’(W) -P |§. + P Re D ’|f.(V^) 

- pR *££S£-l» 4 *“° < 2 - 3 > 

f£ + KV,TVRa|f: + PReU-f£ 

-IV2T* = 0 (2.4) 

where the tilde denotes a dimensionless Laplacian operator, and K denotes a dimensionless 
Jacobian operator. 
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The dimensionless forms lead to two forcing parameters in (2.3) - (2.4). The Reynolds 
number Re is given by 

= ly ( Zr J') I Zrj-i/TTV (2-5) 

and represents dynamic forcing imparted by the basic wind. The Rayleigh number Ra is 
given by 

Ra = Rai + Raf = [(7e _ 7d) z <j> (T 15 — Tsa)] gZipS/ir^TsaWt (2-6) 

and represents thermodynamic forcing. Here, we have separated the thermal forcing into two 
terms in order to define the energetics (section 3) and to aid eventual application of the model 
to the atmosphere: Rai < 0 is the slightly stable boundary layer contribution that is 
proportional to the positive potential temperature gradient or equivalently to the difference 
7d~7e between the dry and environmental adiabatic lapse rates, and Raf is the thermal 
forcing contribution that is proportional to the difference Tib— Tsa between the lower 
boundary and surface air temperatures. 

Three other dimensionless variables appear in the system (2.3) - (2.4). The eddy 
Prandtl number P = ujn is the ratio of the constant eddy viscosity v and the constant eddy 
thermometric conductivity /c; as noted by Laufersweiler and Shirer (1989), this assumption is 
reasonable for modeling boundary layer rolls. Although difficult to estimate, we may use the 
results of Brummer (1§85) to guide our choice of a value for v\ normally atmospheric values of 
P are assumed to be near 1. The roll aspect ratio a is defined as 

2Zrp 

a = — — (2.7) 

L 

in which the domain height is generally chosen to be the cloud top or inversion base, and 
the roll wavelength L is obtained from the value of a. The variable U (z*) represents the 
dimensionless cross-roll wind profile. Owing to the dimensionless forms chosen, we have the 
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constraint | U*( 7 r) | = 1. In this system, the roll-parallel wind component V (z*) has been 
decoupled via the elimination of the Coriolis parameter and the assumption of two- 
dimensionality. Thus we are not considering some longitudinal contributions that may be 
important. 

Using the formula, 

U(z) = U s (z) sin(/3) - V s (z) cos(0) (2.8) 

we may rotate the cross-roll winds into a standard coordinate system for which the eastward 
U s (z) and northward V s (z) wind components are labeled with the subscript s (Fig. 1). Here 0 
is the roll orientation angle that is defined to be the angle between a standard reference 
direction x s and the roll axis y; thus, for example, when 0 = 0* , U = -V s . Positive values are 
measured counterclockwise, negative values clockwise. For a given wind profile, the values of 
a and 0 represent geometric characteristics of the modeled roll circulations. Because the 
cross-roll wind U(z) changes as the value of 0 is changed, the two-dimensional solutions can 
react to the complete horizontal wind profile by choosing an optimal angle 0 P . In this way we 
are able to incorporate some aspects of both horizontal dimensions into the modeled dynamic 
forcing. 

Fourier expansions, composed of temporally dependent amplitudes and spatially 
dependent trigonometric functions, are used to represent the dependent variables ip* and T : 

tf*(x*,z*,t*) = #(x*,z*,t*) + ^(z*,t*) 

= [^i(t*) sin(x*)sin(qz*) + ^(t*) cos(x*)sin(qz*) 

+ V» 3 (t*) sin(x*)sin(nz*) + ^ 4 (t*) cos(x*)sin(nz*)] 

+ [^s(t’) sin((q-n)z») + ^ 6 (t*) sin((q+n)z*)] 


(2.9) 
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T*(x*,z*,t*) = T*(x*,z*,t*) + T b (z*,t*) 

= [Tj(t*) cos(x*)sin(qz*) + T 2 (t*) sin(x*)sin(qz*) 

+ T s (t*) cos(x*)sin(nz*) + T^t*) sin(x*)sin(nz*)] 

+ [T 5 (t*) sin((q-n)z*) + T 6 (t*) sin((q+n)z*) 

+ T 7 (t*) sin(2qz*) + Tg(t*) sin(2nz*)] (2-10) 

The roll solutions ^ and T* are given by the eight terms involving spectral components V'i 
through ipi and Ti through T4, while the nonlinear modifications V’fc and T b of the initial 
wind and temperature profiles are given by the six horizontally constant terms involving 
components ^5, tpf, and T5 through T 8. To permit horizontal roll propagation, both sine and 
cosine functions of x* are required (Pyle 1987), while to satisfy the vertical boundary 
conditions, only sine functions of z* are used. Two general wavenumbers, q and n, are 
included in the vertical representation of (2.9) - (2.10) for study of the inflection point 
instability and for improved representation of the thermal instability, which requires at least 
one harmonic (Stensrud 1987a). Although the above truncations are sufficient for the 
approximation of most simple flow patterns, more spectral modes would be needed to 
quantify completely the initial profile U(z) and the solutions at larger supercritical values of 
the forcing. In section 5, we verify that the Ekman profile we study produces stability results 
that are not sensitive to increased vertical resolution in the model. 

Upon substituting the expansions for ^ and T into (2.3) — (2.4) and integrating 
appropriately over the domain, we obtain the 14 time-dependent spectral equations given in 
Appendix A. We note that the nonlinear terms in the spectral equations (A5) - (A6) and 
(All) - (A14) for the profile modification coefficients correspond to the usual relations 

= -du*w*/dz* and cfTl/dt* = -dw*T*/dz* } where the overbars denote horizontal 
averages (e.g. Chlond 1987). Definitions for each of the coefficients ai, bi, Ci and di in these 
equations are shown in Table Al. Coefficients ci, multiplying the dynamic forcing parameter 
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Re, contain the Fourier coefficients Aj and Tj of the dimensionless cross-roll wind profile 
U*(z*). The Fourier integrals are defined in Table A2. Values for Ai and Ti may be obtained 
from either idealized or observed wind profiles. Here, idealized results using the Ekman wind 
profile are analyzed, while in a companion paper, observational results from the 1987 
stratocumulus experiment FIRE are considered; a preliminary report of these results is given 
in Shirer and Haack (1990). 

3. Energetics 

The dimensionless Boussinesq system (2.3) - (2.4) contains sources and sinks of both 
available potential AE and kinetic KE energies that contribute to the growth and 
development of roll circulations. Sources of roll energy are of both thermal and dynamic 
type, while sinks include eddy dissipation and roll modifications of the initial background 
profiles. These individual energy contributions and the interactions between the secondary 
and background flow may be analyzed by separating the dimensionless energetics 
components into four parts: roll kinetic KER and available potential AER energies, and 
background kinetic KEB and available potential AEB energies. Fundamentally, the 
available energy definitions change with the sign of Ra = Rai + Raf in (2.6). 

A co mm on environment for secondary instability is produced when cold air overspreads 
warmer water. These conditions create strong capping temperature inversions for which a 
constant domain height z^, is appropriate. Typically in this case, the boundary layer has a 
slightly stable initial lapse rate (Rai < 0), and buoyant forcing at the sea surface (Raf > 0). 
When | Raf | > | Rai | , the value of the total thermal forcing rate Ra is positive, and 
contributions from both thermal and dynamic sources can force secondary instabilities; 
conversely, when | Raf | < | Rai | , the value of Ra is negative, and only dynamic sources can 
force the instabilities. 
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Here we examine the case Ra > 0 for which both roll thermal generation and AER/KER 
conversion, or heat flux, terms are obtained. Although similar formulations of the energetics 
have appeared elsewhere (e.g. Kaylor and Faller 1972; LeMone 1973), we include our version 
here to make clear the sources and sinks of energy in the present model, with particular 
attention given to the origins of the profile modification terms. The appropriate definitions 
for KER, AER, KEB and AEB are 


KER = j j 27r J T |V^| 2 dx*dz* 


AER = 


if 2 */* P T* 2 dx*dz* 
^ J 0 J 0 1 


FR5) 

KEB=ij?V [RePU‘-|$] J dx>dz 


0 ' 0 




(3.1) 

(3.2) 

(3.3) 

(3.4) 


in which -difc/d z* = ug is the modification to the initial wind profile by the secondary flow, 
and Tf is the dimensionless form for Tf in (2.2). The definitions for KEB and AEB are valid 
provided that | difo/ dz*\ < Re P | U* | and | T^ | < | Tf | , which we confirmed for the cases 
examined in sections 5 and 6. The following Ra > 0 energy rate equations are obtained from 
the partial differential system (2.3) - (2.4): 


(HF) (RS) (K r -MOD) (K r -DIS) 

KER = P *|g + Re|& ||j f£ -p# K“( $ 5 . | V 2 1 ’, | ’] dx-dz- 


(3.5) 


-(HP) 


(GA) 


{A r -MOD) 


aer = p jf j^[-t + .Ml T -=pi; T r K ” < 16 ; ' T b) - 


(A r - DIS) 


dx*dz* 


(3.6) 
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-(RS) 


(K b -MOD) 


(Kb-DIS) 


KEB = P f 27r J 7r 
} 0 0 


~ Re l|s I v 2 1 2 dx*dz* 


(3.7) 


-(GA) 


(Ab-MOD) 


(Ab-DIS) 


AEB=p/f^[-_g;T;^-^T; iKr 


dx*dz* 


(3.8) 

When Ra < 0, we redefine AER and AEB in a typical way so that (-Ra) replaces (-Rai) in 
the denominators of (3.2) and (3.4). In this case, the available potential energy rate 
equations simplify so that terms -(HF) and (GA) combine to form only a heat flux term 
-(HF) in (3.6), and the term -(GA) is eliminated from (3.8). 

In Appendix B, the energy rate equations (Bl) - (B6) corresponding to the spectral 
system (Al) - (A14) are given for the case Ra > 0. These rate equations contain spectral 
representations of the terms labeled in (3.5) - (3.8) for the partial differential system. Since 
we retain only a few Fourier coefficients in the variable expansions (2.9) - (2.10), the 
representations of these energy terms are only accurate to first order. As in (3.1) - (3.4), 
energy contributions generated by the roll perturbations are separated from those owing to 
the background flow so as to elucidate the nonlinear modifications of the initial 
environmental wind and temperature profiles. Each term contributes to the complete energy 
budget of the system, as is shown schematically in Fig. 2. In this diagram, the four energy 
pools are depicted, and in Table 1, each of the energy sources and sinks is described. 

From Fig. 2 and Table 1, we see that rolls may extract background energy in the usual 
mann er via both mechanical generation (RS) and, in statically unstable boundary layers, via 
thermal generation (GA); roll available potential energy (AER) is converted to roll kinetic 
energy (KER) in both stable and unstable boundary layers via vertical heat flux (HF). 
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Conversely, the secondary flow may alter the initial background state via terms (Kb-MOD) 
= -(Kr-MOD) and (Ab-MOD) = -(Ar-MOD); typically, these terms are not displayed 
because the roll and roll modification components of the response are combined. Finally, 
terms (K-DIS) and ( A-DIS) represent dissipation from each of the four pools of energy. The 
rate of change of total roll energy E corresponding to the spectral system (Al) - (A14) may 
be written as 

• • • 

E = KER + AER 

= [(RS) + (GA)] + [(Kf-MOD) + (Ar-MOD) + (Kr-DIS) + (Ar-DIS)J (3.9) 

For steady energetics, bifurcations to temporally periodic roll solutions occur when the 
energy sources, given by (RS) and (GA), balance the energy sinks, given by (Kr-MOD), 
(Ar-MOD), (Kr-DIS) and (Ar-DIS). Values for E and for each of the individual energy 
terms are calculated in section 5b2 from solutions to the nonlinear spectral model. 

In the following section, we outline how we locate in (Ra,Re)-parameter space the 
transition curves representing each of the possible roll instability modes. We also discuss the 
application of a higher resolution linear model used in the verification of the two-vertical 
wavenumber spectral model results. 

4. Linear stability analyses 

In this section we first outline the analytical linear stability analysis that is used to 
approximate the onset of roll modes in our spectral model, and then we discuss the more 
accurate numerical stability analysis of the linearized Boussinesq equations that is used to 
find an Ekman profile whose stability properties can be captured relatively well by our 
model. As noted by Stensrud (1987a), in both cases the roll modes of interest are found by a 
linear analysis of the motionless conductive solution. 
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a. Bifurcation equation for the roll modes 

As we noted in the introduction, a major advantage of use of low-order spectral models 
is that many of their solutions can be obtained analytically, thereby allowing a thorough 
exploration of roll behavior in forcing parameter space. Although linear stability analysis is a 
commonly used technique, the determination of the bifurcation equation governing the onset 
of the roll modes is somewhat complicated, and so we review briefly how we find this 
equation. 

The linear system obtained from the spectral model (Al) - (A14) contains six equations, 
originating from (A5) — (A6) and (All) — (A14), that are purely dissipative and so do not 
contribute to a change in stability. The remaining eight may be written as two sets of four 
complex equations in which one set is the conjugate of the other. Thus, only one set of four 
equations need be analyzed since it contains the same stability information as the other set 
(Pyle 1987). The form exp[At*] is assumed for a perturbation expressed in the new complex 
variables ip j and Tj; the real parts of the eigenvalues A indicate growth or decay of this 
perturbation and the imaginary parts represent its temporal periodicity. A nontrivial 
solution is found when the determinant of the matrix of coefficients vanishes, which in this 
case produces a complex fourth-order characteristic equation in A (Haack— Hirschberg 1988). 
Transitions to temporally periodic solutions are found by assuming that the real part of A 
vanishes and that the imaginary part of A is the dimensionless limiting frequency u* of the 
branching solution. These Hopf bifurcation points (e.g. Marsden and McCracken 1976; Pyle 
1987) indicate critical values Ra c and Re c of the forcing rates beyond which propagating roll 
solutions develop. Minimum values of Ra c and Re c provide the smallest forcing values 
needed for secondary instability and are presumably the first values reached as the 
atmosphere becomes unstable. 
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Substitution of A = iw* into the complex fourth-order characteristic equation yields two 
real equations in w*, Ra c and Re c . Upon combining these equations via the method of 
eliminants (Stensrud 1987b), we remove the frequency and obtain a Hopf bifurcation 
equation that is of sixth-order in Ra c (Haack-Hirschberg 1988). The coefficients of this 
equation were determined analytically with the program FORMAC on the PSU IBM 
3090-600s. They are lengthy functions of Re c , aspect ratio a, orientation angle p, vertical 
wavenumbers q and n, Prandtl number P, and Fourier coefficients Aj and T i of the basic 
wind. For a given arbitrary wind profile, we choose values for q, n and P and cycle over 
appropriate ranges for Re c , a and p. Minimum values of the critical forcing rates Ra c and Re c 
are found and these give preferred roll alignments P p and aspect ratios a p . From values of 
(Rac)min, (Re c )min, Pp and a p , preferred values of the (real) dimensionless frequency 
magnitude | u* \ are also calculated; the magnitudes of the limiting dimensionless 
propagation rates of the rolls are given by | c* | = | lj* | . Using the dimensional form u = 
[(2t 2 /c)/(z I j,L)]o^ of the frequency, where /c is estimated to be 25 m s /s (e.g. Shirer and 
Brummer 1986), and the definition | c | = (L/2x) | a/| , we may obtain a preferred value | c p | 

= (7rAc/z,p) | w* | for the phase speed. When these Hopf bifurcation points are displayed in 
(Ra,Re)-parameter space, transition curves for the various roll modes are produced (e.g. 
Shirer 1986; Stensrud 1987a); in section 5 we find that three principal modes are possible. 

b. Numerical analysis of the Boussinesq equations 

Although the bifurcation equation discussed above yields analytic approximations of the 
critical forcing rates and the preferred response parameters for the various roll modes, it only 
incorporates the effects of five of the Fourier coefficients of the initial background wind 
profile (cf. (2.1)). Because information in other Fourier terms of this wind profile may 
significantly affect the values of the forcing rates or the response parameters, a more 
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thorough stability analysis must be done to ensure that the two-vertical wavenumber model 
(Al) - (A14) may be used successfully to examine the possible secondary flows created by a 
particular wind profile. This more detailed analysis can not be done analytically, thereby 
making a complete survey of parameter space tedious and computationally expensive. 

The stability results produced by spectral expansions having two (2— WN), three 
(3-WN), and four (4-WN) vertical wavenumbers are obtained numerically and compared 
with those given by the Hopf bifurcation equation for the two-wavenumber model. The 
Fourier expansions of the basic wind profile have ten terms in the (3-WN) case and 17 terms 
in the (4— WN) case. Many of the terms in the 8 x 8, 12 x 12 and the 16 x 16 stability 
matrices, which are associated with the (2— WN), (3-WN) and (4— WN) representations 
respectively, contain spectral coefficients of the initial wind profile. However, because these 
coefficients involve at most only two of the vertical wavenumbers q and n (see Table A2), it is 
easy to specify the terms in these determinants merely by considering all combinations of the 
relevant wavenumber pairs in the general forms given in Haack-Hirschberg (1988). Because 
we normally find that the smallest wavenumbers give the best representation of an initial 
wind profile, we limit our more detailed analyses to combinations of q = 1 and n = 2 in the 
(2-WN) case, q = 1,2 and n = 2,3 in the (3-WN) case, and q = 1,2,3 and n = 2,3,4 in the 
(4-WN) case, where q # n. 

To find the approximations of the roll mode transition curves given by the (2-WN) and 
the higher resolution (3-WN) and (4-WN) representations, we must find all the eigenvalues 
of the above matrices for a large number of values for the parameters Ra, Re, a and /?. 
Specifically, a point (Ra,Re) on the transition curve as well as the associated preferred values 
of a, (3 and u* are obtained as follows: The magnitude of Re is fixed and the value of Ra is 
gradually increased. For a particular value of Ra, the values of a and /? are varied and the 
eigenvalue having the largest real part is saved in each case. The pair (a,/?) is found for which 
this largest real part reaches a maximum value. A preferred bifurcation point on the 
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transition curve is obtained once the magnitude of Ra is increased sufficiently that this 
maximum value is zero. 

If the values of the critical forcing and preferred response parameters do not vary 
significantly as additional wavenumbers are added to the problem, then we may conclude 
that a good representation for the instability modes has been obtained with only two vertical 
wavenumbers. Although to validate a particular wind profile we must investigate a large 
number of parameter values, we are guided somewhat in our choice of parameter ranges by 
the analytic results given by the Hopf bifurcation equation for the two-wavenumber model. 

5. Model results using the Ekman profile 

To facilitate comparisons with previous studies of idealized flow, we consider the 
stability results and nonlinear solutions forced by an Ekman spiral having Ekman depth D, 
where xD /4 is the height of maximum shear in the boundary layer. We define the 
dimensionless eastward and northward wind components for this profile as (Shirer 1986) 

u s( z *) = I Y*| [ 1 -exp(-z*/D*) cos(z*/D*)] (5.1) 

V I( Z *) = I Ygl exp(-z*/D*) S in(z*/D*) (5.2) 

where D* = Di/z^, and where | V* | is the magnitude of the dimensionless geostrophic wind 

£ £ * 
given by the constraint U s (7r) 2 + V s (tt) 2 = 1. We show below that an Ekman depth near D 

= 1 produces consistent stability results when more wavenumbers are included in the 

analysis. To obtain D* = 1, we choose D = 190 m and Zrp = 600 m; such a value of z,j, is 

typical of marine boundary layer depths (e.g. Brummer 1985). Because the Ekman depth D 

is well below the inversion height z^,, use of constant eddy viscosity v and conductivity k is 

appropriate. We note that other authors (e.g. Faller and Kaylor 1966, 1967; Brown 1972; 
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Asai and Nakasuji 1973) used infinitely deep domains and determined roll circulation depths 

* 

from the solution itself. Despite these differences, we find below that our results using a D = 
1 Ekman profile agree rather well with these published results as well as with observations 
reported by LeMone (1973), LeMone and Pennell (1976) and Brummer (1985). 

Using the definitions in Table A2 with the above noted values q = 1 and n = 2, we 
calculate the values of the spectral coefficients Ai and Ti of the Ekman wind (5.1) - (5.2); we 
also calculate the coefficients for the wavenumber pairs (q,n), q ^ n, given by q = 1,2,3 and n 
= 2,3,4 that are needed in the validating numerical stability analysis. When the 
dimensionless form of (2.1) is substituted into the Fourier coefficient definitions (Table A2), 
expressions for the Uj may be obtained to produce the following (2— WN) representation: 

U (z*) = (Ai H- ^ r i + r 3) — f 1 3 sin(z*) + ^{4 A2 + Tj) cos(z*) 

-jj(A 2 + r 2 ) sin(2z*) + ^{r, -|r,)cos(2z-) (5.3) 

We show in Fig. 3 the Fourier approximations calculated from (5.3) for the along-roll V = 
U* and cross-roll U* = -V* wind components when /?= 0* (cf. Fig. 1). The solid curves 
represent the original Ekman wind profiles given by (5.1), (5.2) and (2.8), and the dashed 
curves show the approximated profiles given by (5.3). These results indicate that an 
excellent representation of the Ekman wind can be obtained with a Fourier series having as 
few as two vertical wavenumbers and five Fourier coefficients. 

o. Transition curves 

Having calculated values of the Fourier coefficients that reproduce the Ekman wind 
profile, we next use the methods described in section 4 to perform linear stability analyses of 
the conductive solution. Here and in the remaining sections, we set typical average values of 
u = k = 25 m 2 /s (P = 1) based on computations by Shirer and Brummer (1986). To obtain 
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the minimum values of the forcing rates, we first find the minimum values of Ra c and Re c 
given by solving the Hopf bifurcation equation of section 4a for appropriate ranges of the 
other parameters. The values of the aspect ratio a and the orientation angle p are varied first 
using a coarse resolution, from 0.1 < a < 2.0 in increments of 0.1 and from —90 < P < 90 in 
increments of 10* , and then on a finer grid in which a is incremented every 0.05 and P every 
degree. These results are then used to guide the higher resolution numerical analysis of 
section 4b. 

For the case D* = 1, we compare in Fig. 4a the transition curves for the two- 
wavenumber (2-WN) (solid), three-wavenumber (3-WN) (dashed) and four-wavenumber 
(4-WN) (dotted) approximations obtained using the numerical method described in section 
4b. It is this value of Ekman depth that yields the closest agreement among the three 
transition curves. Along these curves, the preferred values of a, P, and u* are given. We see 
dearly from the figure that the transition to secondary flow is well represented by the 
(2-WN) curve when Re is less than 60 or so, but that this curve begins to depart rapidly from 
the other two as the magnitude of Re is increased above 60. The (3-WN) and (4-WN) curves 
remain dose for wider ranges of Ra and Re. This result was found generally for a large 
number of Ekman depths and implies that a (3-WN) modd may be more generally applicable 
(e.g. Stensrud and Shirer 1988). In the range 0 < Re < 60, the orientation angle results among 
the three approximations differ by at most 4* and the aspect ratio values by no more than 
0.1, indicating very consistent results. Moreover, the response parameter values themselves 
compare well with those found by others who theoretically studied Ekman flow (e.g. Lilly 
1966; Faller and Kaylor 1966; Brown 1970, 1972; Etling 1971; Asai and Nakasuji 1973), and 
by those who observed the atmospheric boundary layer (e.g. LeMone 1973; Walter and 
Overland 1984; Briimmer 1985). For example, for modes near neutral stratification (Ra ~ 0), 
we obtain the values a p ~ 0.6, L p ~ 2000 m ~ 10D, P» 10’ and |c p | ~ 0.2 m/s that are 
consistent with those reported by the above authors. 
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We conclude that for this Ekman profile, we may use the model (Al) - (A14) to study 
the modes developing when Re c < 60. However, we emphasize that this cutoff value is not a 
general one; each wind profile must be checked individually. We may also conclude that if 
the preferred orientation angles produced by using a roll-modified Ekman wind profile as the 
basic state differ by more than 4* from those given by the pre-roll profile, then we have 
evidence that the errors obtained using observed profiles as basic states might be more likely 
attributable to the contamination of the observations by the roll solutions than to the limited 
truncation used to develop the model (Al) - (A14). We return to this point in section 6. 

In Fig. 4b we show the transition curves obtained from the bifurcation equation for the 
two-wavenumber spectral model, as described in section 4a. Aside from minor differences in 
frequency — which are more accurately given by the bifurcation equation — these values are 
equivalent to those given in Fig. 4a for the (2-WN) numerical analysis of the linear model. In 
Fig. 4b we separate the transition curve into the two principal modes: the thermal-q mode 
(dashed curve) that produces a solution dominated by wavenumber q = 1 and is associated 
with Rayleigh-Benard convection when U(z) = 0, and the inflection point mode (solid curve) 
whose transition curve passes through the neutral stratification value Ra = 0. A third mode 
was also found at large values of the Reynolds number (Re > 100) and in strongly stable 
stratification (Ra < 0); however, comparison of this dynamic mode with ones produced by 
the higher resolution models revealed that, while appearing in all three analyses, it 
nevertheless can not be well represented by our model. Thus, we restrict attention to values 
of Ra near those for neutral and unstable stratification. 

As in Fig. 4a, the preferred values of the response parameters are given. We note that 
with increasing values of Re c , the thermal— q mode is replaced by the inflection point mode 
near neutral stratification (Ra = 0). For a variety of wind profiles, a smooth transition from 
the thermal-q to the inflection point mode occurs, indicating a possible link between these 
two instability mechanisms. Other thermal modes associated with wavenumbers n = 2 and 
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combinations of q and n are also common, but are not shown here because they occur at 
larger values of the forcing parameters and are therefore of little interest. Hereafter, we refer 
to the thermal-q mode as simply the thermal mode. 

Two paths denoted by arrows and labeled with points A and Bi are drawn in Fig. 4b; 
these paths, which we consider below, represent two possible evolutions of the atmosphere as 
the values of Ra and Re increase. In each case, the value of Ra > 0 is a statically unstable one 
so that the resulting circulation is both thermally and dynamically forced, and the value of 
Re < 60 is in the range for which the model is valid (Fig. 4a). We summarize in Table 2 the 
minimum values (Ra c ,Rec)niin of the critical forcing rates needed for roll development and the 
supercritical values (Ra c ,Re c )sup of the forcing rates chosen at the labeled points. In section 
6, we ascertain how the rolls modify the initial profiles of wind and temperature by examining 
the effects of increasing the buoyant forcing, as given by increasing the value of Ra 
incrementally from point Bi to B2 to B3 in Fig. 4b. In the following subsection, the solutions 
corresponding to points A and Bi are used to compare the roll stream function and 
perturbation temperature patterns, the energetics profiles, and the vertical fluxes of heat and 
momentum produced by the inflection point and thermal instability modes. 

b. RoU. solutions 

To e xami ne the nonlinear solutions given by the values of the forcing rates at points A 
and Bj in Fig. 4b, we choose values of a and /? near the preferred values on the transition 
curve closest to each supercritical point (see Table 2). We then numerically integrate the 
nonlinear equations (Al) - (A14). The rate of change of total roll energy is calculated from 
(3.9) using the definitions in Appendix B, and in all cases we find that energetically steady, 
temporally periodic roll solutions and steady background modifications occur. As noted 
earlier, the decay of the perturbations to a stable solution typically requires less than four 
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hours, which justifies the neglect of the Coriolis terms. Figure 5 shows a schematic 
representation of the stream function pattern produced for a small, positive value of the 
preferred orientation angle /? p . In this case, the roll circulations tilt downstream and 
propagate at a fixed amplitude in the -x direction. The basic structure and tilt of these rolls 
are consistent with those shown, for example, by Faller and Kaylor (1966), Brown (1970, 
1972) and LeMone and Pennell (1976). 

1). Stream function and perturbation temperature patterns 

At a fixed time, the dimensionless roll stream function patterns ip* (from (2.9)) that 
correspond to the four points A and Bi are shown in Fig. 6. All patterns propagate from right 
to left in the direction of the cross— roll wind (Fig. 3b). As the dynamic and thermodyn ami c 
forcing rates Re and Ra are increased from point Bi to A (Fig. 6b to 6a), we note that the 
circulation tilt increases downwind in response to the greater wind speeds at the top of the 
domain (cf. (2.5)). The effects on the ip* fields of increasing only the thermal forcing rate Ra 
are shown in Figs. 6b to 6d. Although the circulation patterns remain relatively unchanged, 
the corresponding dimensional maximum upward velocities increase from about 0.1 to 
1.1 m/s. Qualitatively, these rolls have characteristics similar to those given by Brown’s 
(1970) model of neutrally buoyant secondary flow (cf. his Fig. 11, in which propagation 
occurs from left to right) and by Brown’s (1972) model of stratified flow (cf. his Figs. 8 and 9; 
propagation is also from left to right). 

The dimensionless roll perturbation temperature patterns T* (from (2.10)) that 
correspond to the points A and Bi are shown in Fig. 7. In each case, positive temperature 
perturbations correspond to regions of upward motion in Fig 6. Direct thermal circulations 
are expected for cases of statically unstable stratification and produce the positive vertical 
heat fluxes shown below in Fig. 10a. As the thermal forcing rate Ra is increased from point 
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Bi to B3 (Fig. 7b to 7d), more efficient upward transports of heat are obtained by the reversal 
in the tilt of the Tr field from upwind to downwind. Owing to the Boussinesq approximation, 
the temperature perturbation in our model best corresponds to either the density 
perturbation or the potential temperature perturbation in observations; the tilt in the 
temperature perturbation patterns in case B3 most closely matches the observed cross 
sections through thermally forced rolls given by LeMone and Pennell (1976; cf. p v in their 
Fig. 12) and Briimmer (1985; cf. 6^ in his Fig. 7b). 

Comparison of the roll solutions with LeMone’s (1973, her Table 2) roll observations 
indicates that the circulation characteristics agree well with her values. LeMone (1973) 
estimated the maximum roll circulation speeds Umax to be between 0.077 1 V g | and 0.15 1 V g | . 
Using appropriate values in the expressions (2.5), (5.1) and (5.2), we obtain magnitudes of V g 
for Re = 30 and 60. From the roll solutions we calculate values of Umax = | -d$/ | max to be 

0.07 1 V g | for case A and 0. 19 1 V g | for case B v Moreover, the cross-roll propagation velocity 
c p , which ranges from -0.2 to -0.3 m/s (Table 2), is of the same sign and order of magnitude 
as the boundary layer average cross-roll wind velocity U, in agreement with results usually 
found in observations (e.g. LeMone 1973) and models (e.g. Becker 1987; Chlond 1987). Thus, 
we conclude that the 14-coefficient spectral model represents rather well the essential 
characteristics of near-neutral two-dimensional secondary flow. 

2 ). Energy budgets and vertical transports 

The magnitude and sign of each energy term in (Bl) — (B6) provide useful information 
about the relative importance of the mechanisms responsible for producing roll circulations 
and for modifying the initial wind and temperature profiles. These terms may be calculated 
from the solutions to the integrated spectral equations (Al) - (A14). First, values for the 
thermal surface forcing Raf and the environmental stratification Rai must be obtained from 
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the value of Ra. Using the definition (2.6) and assuming T sa = 16* C and near-neutral 
conditions ( 7 d - 7e = 0.1 * C/600 m), we calculate a value of Rai ~ -10. Then Raf (and so 
Tib) is determined by subtracting Rai from the given value of Ra. 

For points A and Bi in Fig. 4b, we form energy budget diagrams in Fig. 8 showing the 
magnitudes of the sources, sinks and conversions of energy. Both a mechanical source of 
energy, via the Reynolds stress term (RS), and a thermal source of energy, via the generation 
term (GA), contribute to roll development. Consistent with the fact that Ra > 0, in all cases 
the heat flux term (HF) is positive so that energy is converted from roll available potential to 
kinetic forms. Comparing the magnitudes of the energy sources for points A and B i reveals 
that mechanical generation is the larger source term for rolls driven by the inflection point 
instability mechanism (Fig. 8a), while thermal generation is the larger term for rolls excited 
by thermal instability (Figs. 8b,c,d). As the thermal forcing rate is increased from point B i 
to B 3 (Figs. 8b to 8d), we note a marked increase in the magnitude of the (GA) term, as well 
as greater channeling of roll energy to the background flow via the modification terms 
(Kb-MOD) and ( Ab-MOD). These roll modification terms become important at larger 
supercritical forcing rates. In case B 2 for example, the ratio of the modification terms to the 
source terms is roughly one— half. Thus we might expect that significant modification to the 
initial wind and temperature profiles would be possible once the thermal forcing rate is 
increased sufficiently. 

Figure 9 shows vertical profiles for every term contributing to the energy budgets in 
Figs. 8a,b,c. Solutions for points A (inflection point mode), Bj and B 2 (thermal modes) 
correspond to dashed, solid and dotted line types respectively. Energetics profiles for the 
supercritical point B 3 are omitted since they contain shapes consistent with those shown for 
B! and B 2 , but with larger magnitudes. Asai and Nakasuji (1973) present energetics profiles 
in their study of the stability of the Ekman boundary layer. Although their upper boundary 
is infinite, vertical profiles of Reynolds stress, generation and dissipation terms, 
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corresponding to the inflection point and convective modes, agree well with those shown here. 
Brummer’s (1985) energetics analyses of KonTur data include Reynolds stress (RS) profiles 
that indicate a roll kinetic energy source in the lower half of the domain near the altitude of 
the inflection point in the cross-roll winds, as well as generation (GA) profiles that reveal roll 
available potential energy in the middle of the domain where the magnitudes of the 
perturbations are greatest. Both of these results are in agreement with those shown here in 
Fig 9. 

/ 

The dimensional profiles of the horizontally averaged, roll-induced vertical heat (wT ) 
and momentum (uw) fluxes are shown in Fig. 10. These flux profiles are consistent with 
those given by Asai and Nakasuji’s (1973) study of Ekman flow, by Etling and Raasch’s 
(1987), Becker’s (1987) and Chlond’s (1987) higher resolution numerical studies of boundary 
layer circulations and by LeMone’s (1973) and LeMone and Pe nn e l l’s (1976) observational 
studies of roll vortices. That is, in a statically unstable environment, positive heat fluxes 
represent a stabilizing process in which the roll perturbations transport relatively warmer air 
upward and cooler air downward. Although most investigators find positive cross— roll 
momentum fluxes, our results in fact have a familiar form that is consistent with the 
downwind tilt of the roll circulations (Figs. 5 and 6). Because the changes d\J/dtin the 
initial cross-roll wind are given by dll/dt = -d(uw)/dz (e.g. Chlond 1987), the momentum 
flux profile should lead to an increase in the cross-roll wind speed below 300 m and a decrease 
above this level. This corresponds to down-gradient transports of cross-roll momentum, 
consistent with the results of Faller and Kaylor (1967) and Brown (1970), and with the 
results of section 6. 
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6 . Modifications to the basic state profiles 

We noted earlier that a fundamental characteristic of nonlinear roll circulations is that 
they can significantly alter the basic state profiles. Faller and Kaylor (1967) and Brown 
(1970) produce Ekman spirals that have been modified by the secondary flow. According to 
Brown (1980), roll modifications may change the basic flow by as much as 10 to 20%. In each 
case, the altered flow has maximum velocities and is more closely aligned with the direction 
of the geostrophic wind at low levels, and becomes supergeostrophic at upper levels (Brown 
1970). These results suggest that cross-roll modifications are channeled via the Coriolis force 
into the along-roll wind component, thereby increasing the roll-parallel flow and reducing 
the roll-perpendicular flow. In contrast, the numerical studies of Chlond (1987) and Etling 
and Raasch (1987) indicate that modification occurs predominantly in the cross-roll 
component of the basic wind profiles. Each of the above models differs from the present one 
by the inclusion of the longitudinal velocity component and Coriolis turning. Thus, their roll 
perturbations may alter the roll-parallel component of the initial basic wind profile, while 
only the roll-perpendicular component may be altered in the present model. We are 
therefore examining the roll modifications that can occur on a short time scale, before 
Coriolis turning can become significant. We find below that with increased supercritical 
forcing rates, the secondary flow decreases the cross-roll shear near the height of the 
inflection point. 

For small supercritical forcing rates, such as those given by points A and Bi in Fig. 4b, 
the dimensional modifications to the initial background wind and temperature profiles are 
insignificant since they change these profiles by only a few percent. However, for the points 
B 2 and B 3 that are much farther to the right of the transition curve in Fig. 4b, the spectral 
components V's, V's and T 5 through Tg that represent modifications to the basic state increase 
significantly in magnitude. From points B 1 to B3, the percentage of energy input that is 


28 


converted to background energy, via terms (Kb-MOD) and (Ab-MOD), grows markedly, 
and so the rolls may begin to noticeably change the mean flow (Fig. 8 ). Although for short 
time scales the neglect of Coriolis forcing is a valid approximation, solutions far away from 
the transition curve may require more spectral modes than we have included in order to 
accurately resolve the nonlinear interactions that lead to a modified basic state. Therefore, 
we only examine the general structure and trends in these modified profiles, with emphasis on 
their potential effect on the results of a linear stability analysis. 

Shown in Fig. 11 are the dimensional profiles of the perturbation temperature 
/ 

modification T b (z) (from (2.10)) for the solutions at points Bi, B 2 and B 3 (solid, dotted and 
dashed curves) on the lower path in Fig. 4b; because the initial temperature profile has a 
near— neutral lapse rate of approximately 9.6 * C/km, the temperature modification shown is 
equivalent to that for potential temperature. As expected, the roll perturbations act to 
neutralize the interior portion of the unstably stratified boundary layer by creating more 
stable background temperature profiles; in each case the lower portion of the domain is cooled 
and the upper portion is warmed. Corresponding to case B 3 in Fig. 11 , we use the 
roll-modified background lapse rate near z = z ^/2 to calculate the modification to the 
thermodynamic forcing rate by the rolls. From its original value of Ra = 50, the thermal 
forcing rate is decreased to a near-neutral value of - 1 . This result is consistent with the 
occurrence of direct thermal circulations and positive heat fluxes associated with thermally 
forced rolls. However, the magnitude of this temperature modification is not likely to be 
measurable in the atmosphere nor great enough to affect the linear stability results given in 
section 5a. 

Conversely, alterations to the basic flow by the rolls do appear to be significant. The 
effects on the basic state wind profile of increasing the supercritical forcing rates are shown in 
Fig. 12 . These profiles represent dimensional profiles of the initial background cross-roll 
wind U(z) (from (5.1), ( 5 . 2 ), (2.8)), the wind modification u b = -dif^/dz (from (2.9)), and 
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the modified cross-roll background wind U + u£ for the solutions at points Bi, B 2 and B 3 . 
We recall that the cross-roll component is determined from (2.8), where 0 = 6" is the 
preferred roll orientation on the transition curve near point Bi (Table 2). As shown in Fig. 
12b, the roll perturbations yield profiles that have increasingly more negative slopes as the 
supercriticality increases from case Bi to B 3 ; the mean wind direction is not changed, 
however. When summed with the original cross-roll Ekman profiles given in Fig. 12a, the 
modified background wind profiles in Fig. 12c are produced. As mentioned above, we 
conclude that these profile modifications are consistent with the momentum fluxes given in 
Fig. 10b. For the forcing rates corresponding to case B 3 , the roll-induced shear is six times 
that of case B i and leads to substantial modification of the original basic flow. Consequently, 
the wind structure throughout the domain can be altered markedly by the roll circulations, 
even in the absence of the Coriolis force. 

In Fig. 13, the dimensional original (solid) and case B 2 modified (dotted) Ekman spirals 
in roll coordinates are shown. As found by Faller and Kaylor (1967) and Brown (1970), the 
cross-roll wind speed of the modified profile is reduced significantly, so that the turning angle 
in the modified Ekman spiral is decreased. In addition, the modified spiral deviates less from 
the geostrophic wind at low levels than does the initial wind, producing a result that is also 
consistent with those of Brown (1970, his Fig. 14). 

In theory, the wind profiles in Fig. 12c correspond to roll-modified winds that might be 
measured during cloud street observations. When values for the Fourier coefficients are 
calculated from these profiles, considerably different stability results might be obtained than 
those shown in Fig. 4b, for which the original, pre— roll Ekman profile was used. We e xami ne 
this possibility by using the roll-modified background wind profile associated with point B 3 
to obtain the transition curves and preferred values of roll characteristics a p , 0 P , v* that are 
shown in Fig. 14. Comparing Figs. 4b and 14, we note that only a thermal mode (dashed 
curve) may excite roll development in the range 0 < Re < 60 for which the model is valid. 
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Owing to the decreased shear near the inflection point in the modified Ekman profile, the 
inflection point mode disappears completely from the stability diagram. More importantly, 
the values of preferred orientation angle /? p differ between 10* and 14* from those shown in 
Fig. 4b and indicate a larger than expected deviation (18* ) from the geostrophic wind 
direction. Such large changes are not entirely related to the use of a limited truncation 
because they are somewhat bigger than the changes shown in Fig. 4a between the (2-WN) 
and (4-WN) stability analyses of the original Ekman profile. In some cases involving 
observed wind profiles taken during cloud street events, use of the roll-modified winds leads 
to errors up to 35* in the estimates of (3 P (Shirer and Haack 1990). Significantly, these 
findings indicate that an appropriate amount of hypothesized roll-induced shear must be 
removed from an observed, cross— roll wind profile in order to determine the initial, basic 
state winds. 

7. Summary and rondusions 

In this study, a 14-coefficient nonlinear spectral model was used to investigate boundary 
layer roll circulations forced by Ekman flow. The Fourier representation of the wind, 
including two vertical harmonics and five Fourier coefficients, approximated the Ekman 
spiral rather well. Decay of small perturbations to stable roll solutions required less than four 
hours, which justified the neglect of Coriolis terms in the model. Thus, longitudinal sources 
of energy, and consequently the parallel instability mechanism, were eliminated. From the 
results of a linear stability analysis, characteristics of the convective and/or dynamic 
instability mechanisms were studied. This analysis was shown to be valid for 0 < Re < 60 by 
comparing its results with those obtained from a numerical stability analysis of linear models 
having three and four vertical wavenumbers in the variable expansions. 

The linear analysis produced transition curves in the forcing parameter space (Ra,Ite) 
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corresponding to each instability mechanism. Here Ra is the Rayleigh number that 
represents thermodynamic forcing from surface heating, and Re is the Reynolds number that 
represents dynamic forcing from the cross-roll wind shear. These curves corresponded to the 
minimum critical forcing rates beyond which stable roll solutions occur. In addition, 
preferred aspect ratios, frequencies, and roll axis alignments relative to the geostrophic wind 
were determined. 

From the analysis of an Ekman wind profile having dimensionless Ekman depth D*= 1, 
we found three possible roll modes. The thermal or convective mode was only preferred for 
cases of unstable stratification (Ra > 0) and weak dynamic forcing (small values of Re). For 
near— neutral stratification (Ra ~ 0), the thermal mode was replaced by the inflection point 
mode. The preferred roll characteristics agreed well with those obtained from other 
boundary layer studies of Ekman flow (Faller and Kaylor 1966; Brown 1970; Etling 1971; 

Asai and Nakasuji 1973), and observational studies of roll vortices (LeMone 1973; LeMone 
and Pennell 1976; Brummer 1985). 

A new dynamic instability mechanism, the shear mode, occurred only in statically stable 
atmospheric conditions and for larger values of Re. Rapidly propagating rolls, with axes 
aligned at large angles to the mean wind direction and with small horizontal wavelengths, 
were excited by the shear instability. However, we concluded from comparisons with higher 
resolution numerical results that at least three vertical harmonics are required to model this 
new mode, and so it was not investigated further. 

At supercritical values of the forcing rates, temporal integrations of the model equations 
yielded nonlinear solutions for the 14 spectral coefficients. These values were used to 
calculate the roll patterns and vertical transports corresponding to roll circulations driven by 
a particular instability mechanism. For all modes, the rolls tilted and propagated downwind 
at a constant rate and amplitude. All roll-induced fluxes were down-gradient and led to 
near-neutral stratification having reduced cross-roll shear. Energetics analyses indicated 
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that roll circulations excited by the inflection point mode were predominantly driven by 
mechanical generation from Reynolds stress, while those excited by the thermal mode were 
driven predominantly by thermal generation from surface heating. 

The alterations of the wind and temperature profiles were represented by six spectral 
coefficients in the model. For positive thermal forcing rates, the rolls modified the basic state 
temperature by cooling the lower, and warming the upper, portion of the domain. 
Dimensionally, however, this temperature change was not likely to be measurable in the 
atmosphere. 

The rolls were shown to significantly alter the initial wind profile in the sense found by 
Faller and Kaylor (1967) and Brown (1970), but via a mechanism independent of the Coriolis 
force. Use of a roll-modified Ekman wind profile gave markedly different stability results 
that yielded errors of order 10* in the preferred roll orientation angle / ? p . These errors are of 
the size reported by Shirer and Brummer (1986), in which the roll modification to the winds 
was not considered. Even larger errors in /? p are possible, as demonstrated by Shirer and 
Haack (1990) who investigate the influence of the rolls on the observed basic wind by 
studying several hypothesized, pre-roll wind profiles for cases obtained during the 
stratocumulus phase of FIRE. 
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Appendix A 

The Nonlinear Spectral Model 

Here we list the 14 equations composing the nonlinear spectral model. The coefficients 
ai, bi, Ci and di are given in Table Al, and the Fourier coefficients Ai and Ti in Ci are defined 
in Table A2. The spectral system is given by 


= — ai ip4ips + + biTi + Re(ciV '2 + c 2 ^4) — P diV’i 

(Al) 

- a 2 ^e - bjTa - Re(c^i + c 2 ^) - P 

(A2) 

= — a3^2^5 + a4^’2V'6 + b 2 T 3 + Re(C3^4 + C4^ 2 ) — P d 2 tf 3 

(A3) 

^ = a 3 ^i^5 - a4^i^6 - b 2 T4 - Re(c 3 V'3 + 04 ^ 1 ) -P d2^4 

(A4) 

3I* = - ^{^4 - Mz) - P dz1>s 

(AS) 

= - ^(^ 1^4 - i>2h) - p d 4 ^6 

(A6) 

= — ae(^ 5 T 4 - i>zTs) + a 5 (^ 3 T 6 -^ 6 T 4 ) + a 7 ^iT 7 + Ra ^>1 

— Re(c5T 2 + C 6 T 4 ) — diTi 

(A7) 

gYl = -a 6 (^4T5-^5T3) + a 5 (^6T3-V'4T6)-a7^2T7-Ra^2 

+ Re(csTi + C 6 T 3 ) — diT 2 

(AS) 

= - a 6 (^ 5 T 2 - 1>i T s ) + a 5 (fcT e - ^ 2 ) + a 8 ^ 3 T 8 + Ra ^3 

-Re(cgT 2 + C 7 T 4 ) -d 2 T 3 

(A9) 

gyj = — a 6 (^ 2 T 5 — if){T i) + as^eTi — t^Te) — a 8 ^ 4 T 8 — Ra ipt 

+ Re(cgTj + C 7 T 3 ) — d. 2 T 4 

(A10) 

j}y£ - ^ 4 T 2 + ^T 4 - ^iT 3 - toTi) - d 3 T 5 

(AH) 

- ^T 4 + fcTj -ill Tj - *T0 - d<T« 

(A12) 

gjJ = ^T 2 -*T,)-(i s T, 

(A13) 
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g^-^V^-^T^-deTg (A14) 

in which Re is the Reynolds number defined in (2.5), Ra is the Rayleigh number defined in 
(2.6), and P = u/kis the Prandtl number. 
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APPENDIX B 
Spectral Energetics (Ra > 0) 

The following rate equations for the roll and background kinetic (KER), (KEB) and 
available potential (AER), (AEB) energies may be obtained from the spectral system (Al) 

- (A14). The individual terms are described in Table 1, and the corresponding rate equations 
for the partial differential system are given in (3.5) -(3.8). 

(HF ) (RS) 

KER = P(Ti^i-T 2 V'2+T 3 V'3-T4V’4) + Re P[A 2 (n 2 -q 2 )(V’ilMM , 3)] 

(K r -DIS) 

+ (Kr-MOD)-^[(a 2 +q 2 ) 2 (V»i J + tM) + (a 2 +n 2 ) 2 (^ 3 2 + ^ 4 2 )] (Bl) 


AER = 


-(HF) (GA) 

— P(Ti^i— T 2 ^ 2 +T 3 ^ 3 — T 4 ^4) + P(T 1 ^i-T 2 V'2+T 3 ^3-T4^4) 

(Ar-DIS) 

+ (A r -MOD)- =R ^- i [(a 2 +q 2 )(T 1 2 +T 2 2 ) + (a 2 +n 2 )(T 3 2 + T 4 2 )] (B2) 


-(RS) 


(Kb-DIS) 


KEB = -ReP[A 2 (n 2 -q 2 )(^4-^3)] + (K b -MOD)-|(2(n-^)^5 2 + 2(n+q)W] 


(B3) 


-(GA) 

AEB = P(Ti^i-T 2 ^ 2 +T 3 ^ 3 -T 4 ^4) + (Ab-MOD) - 

(Ab-DIS) 

- P — [2(n-q) 2 T 5 2 + 2(n+q) 2 T 6 2 + 8q 2 T r 2 + 8n2T 8 2 ] 
— ita i a 


(B4) 


where the nonlinear modification terms are given by 
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(Kr-MOD) = - (Kb-MOD) = 


a s (^5^i^4 - Mato)(q-n) 2 + - ^e^3)(q+n) 2 (B5) 


(Ar-MOD) = - (Ab-MOD) = 

-Rai fcsC^^Ti — Tg^4T 2 + Tg^jTj-Ts^T^ 

+ a 6 (T 5 ^Ti-T5^4T 2 + T5^iTj-T 5 ^2T 4 ) 

+ a 7(T7^iTi — T7^ 2 T 2 ) - &8(Tg^3T3 — Tg^T*)] (B6) 

The coefficients aj in (B5) — (B6) are given in Table Al. The dynamic and thermodynamic 
forcing parameters Re, Rai and Raf are defined in (2.5) and (2.6). The parameter a is the roll 
aspect ratio defined in (2.7), and P = u/nis the Prandtl number. The Fourier coefficient A 2 
is given in Table A2, and q and n are vertical wavenumbers. 
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FIGURE CAPTIONS 

Fig. 1. Rotation of the standard reference coordinate system (x s ,y s ) into the roll coordinate 
system (x,y). The cloud street orientation /? is the angle between the direction x s and the roll 
axis y (after Shirer 1980). 

Fig. 2. Schematic energy budget diagram composed of roll kinetic (KER), roll available 
potential (AER), background kinetic (KEB), and background available potential (AEB) 
energies for the case Ra > 0. The individual terms are defined in the rate equations 
(Bl) — (B6) and described in Table 1. When Ra < 0, the generation term (GA) is eliminated. 

Fig. 3. Dimensionless Ekman along-roll (V* = Us) (a) and cross-roll (U* = -Vs) (b) wind 
components when (3 = 0* for Ekman depth D* = 1. The solid curves are obtained from (5.1), 
(5.2) and (2.8), and the dashed curves from the Fourier representation (5.3). 

Fig. 4 Curves of minimum critical forcing rates (Ra c , Re c )min for the Ekman profile (5.1) 

- (5.2) when D*= 1. Values of the preferred aspect ratio a p , orientation angle /? p (in degrees) 
and frequency magnitude | oj * | (in parentheses) are labeled along each curve. In (a), the 
results using the two-wavenumber (2-WN) (solid), three-wavenumber (3-WN) (dashed) 

. and four-wavenumber (4-WN) (dotted) numerical stability analyses of section 4b are shown. 
In (b), the two-wavenumber results obtained from the Hopf bifurcation equation of section 
4a are shown, with the dashed curve representing the thermal-q instability mode and the 
solid curve the inflection point instability mode. The arrows indicate two possible evolutions 
of the atmospheric forcing rates Ra and Re from subcritical to supercritical values; solutions 
at the points A and Bi are discussed in sections 5b and 6. The corresponding parameter 
values at each point are given in Table 2. 
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Fig. 5. Schematic depiction of a propagating roll stream function pattern in the presence of an 
Ekman wind profile V for a small, positive orientation angle /? p . The axes (x s ,y s ) represent 
the standard east/north coordinate system, and the axes (x,y) represent the roll coordinate 
system. The parameter z^, specifies the domain height, the roll spacing is given by the 
preferred wavelength L p , and V g is the westerly geostrophic wind vector. Propagation is 
from right to left in the direction of the cross-roll flow. 

Fig. 6. Patterns of dimensionless roll stream function ^ for the points labeled A (a), Bi (b), 
B2 (c) and B3 (d) in Fig. 4b. Propagation is from right to left as denoted by the horizontal 
arrow. Parameter values associated with each point are given in Table 2. 

Fig. 7. Patterns of dimensionless roll perturbation temperature T r for the points labeled 
A (a), B! (b), B2 (c) and B3 (d) in Fig. 4b. Propagation is from right to left as denoted by the 
horizontal arrow. Parameter values associated with each point are given in Table 2. 

Fig. 8. Dimensionless energy budget diagrams corresponding to the points labeled A (a), 

Bi (b), B2 (c) and B3 (d) in Fig. 4b. Each term corresponds to a dimensionless energy rate 
defined in (BI) - (B6) and described in Table 1. The dimensionless roll energy is quantified 
by the value of E, and the boxes represent four separate pools of energy: roll kinetic (KER), 
roll available potential (AER), background kinetic (KEB), and background available 
potential (AEB). Parameter values associated with each point are given in Table 2. 
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Fig. 9. Vertical profiles of each of the dimensionless energy rates defined in (Bl) - (B6) for 
the points labeled A, B! and B 2 in Fig. 4b. Dashed curves correspond to the supercritical 
point for the inflection point mode (point A), solid curves for the thermal mode (point Bj), 
and dotted curves for the thermal mode (point B 2 ). The boxes (KER), (AER), (KEB) and 
(AEB) represent the four pools of energy. Parameter values associated with each point are 
given in Table 2. 

Fig. 10. Vertical profiles of dimensional vertical (a) heat and (b) momentum fluxes for the 
points labeled A, Bj and B 2 in Fig. 4b. Line types are as in Fig. 9. Parameter values 
associated with each point are given in Table 2. 

/ 

Fig. 11. Vertical profiles of the perturbation temperature modification T b for the solutions at 
the points labeled B i (solid), B 2 (dotted) and B3 (dashed) in Fig. 4b. Parameter values 
associated with each point are given in Table 2. 

Fig. 12. Vertical profiles of the initial cross-roll wind U (a), the wind modification u b (b), 
and the modified background cross-roll wind U + u b (c), for the solutions at the points 
labeled B 1 (solid), B 2 (dotted) and B3 (dashed) in Fig. 4b. Parameter values associated with 
each point are given in Table 2. 

Fig. 13. Dimensional original and modified Ekman spirals and the geostrophic wind vector 
Vg in roll coordinates. The solid curve is from the original wind components (5.1) — (5.2) and 
the dotted curve from the modified winds shown in Fig. 12c corresponding to case B 2 . 
Parameter values associated with this point are given in Table 2. 
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Fig. 14 Curves of minimum critical forcing rates (Ra c , Re c )min obtained from the Hopf 
bifurcation equation of section 4a for the modified cross-roll Ekman wind profile shown in 
Fig. 12c corresponding to point B3. The dashed curve represents the thermal-q instability 
mode. Values of preferred aspect ratio a p , orientation angle P p (in degrees), and frequency 
magnitude | u* | (in parentheses) are labeled along the curve. Parameter values associated 
with this point are given in Table 2. 
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Table 1 


Summary of energetics terms in (3.5) - (3.8). 

ORIGINATING TERMS SOURCE 


TERM DEFINITION 

HF vertical heat flux or 
AER/KER conversion 

RS Reynolds stress or 
mechanical generation 

GA thermal generation 
(Ra > 0) ; or zero (Ra < 0) 
with reformulated AER and 
AEB rate equations 

Kb-MOD initial basic state 

wind profile modification 

Ab~M0D initial basic state 
temperature profile 
modification 

K-DIS roll and background 
viscous dissipations 

A-DIS roll and background 
thermal dissipations 


'#T*/cbc* in (2.3) and KER 

dfr/dx* in (2.4) 

advection terms involving KER 

U* in (2.3) 

(Tib - T sa ) term of Ra AER 

in (2.4); 

or none - - - 

Jacobian term in (2.3) KEB 

Jacobian term in (2.4) AEB 


last term in (2.3) 
last term in (2.4) 


SINK 

AER 

KEB 

AEB 


KER 

AER 


KER 

KEB 

AER 

AEB 



Table 2 


Parameter values associated with the path arrows shown in Fig. 4b. The 
values for (Ra c ,Re c )nnn denote the minimum, critical forcing rates needed for 
roll development, and the values for (Ra c ,Re c )sup denote the supercritical 
forcing rates corresponding to each of the points A, Bj , B 2 and B 3 on the 
paths. Here L p , a p , /? p , c p , |w*| and T p are the preferred values of 
horizontal wavelength, aspect ratio, orientation angle, propagation velocity, 
dimensionless frequency magnitude and period. The quantity L p /z^ represents 
the ratio of the horizontal to vertical roll scales, and U is the average 
cross- roll wind velocity in the boundary layer. Other parameter values are: 
7 d - 7e = 0.1 *C/600 m, v - k - 25 m 2 /s, P = 1, Zj - 600 m, D = 1 , q = 1 , 
n = 2 , A! = - 0 . 112 , A 2 = -0.094, A 3 = -0.104, = 0.094, T 2 = 0.259, and 

r 3 = 0.242. 


H 

(Ra c jB.6c ) mi n 

(Ra c }R-6 c )sup 

L p (a p ) 

(M 


ES 

VSm 


lit 

■S3 

El 

H81 
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(5,60) 

(10,60) 

1.8 (0.65) 6* 

- 0.3 

(2.2) 

1.8 

3.1 

- 0.6 

Bj 

(10,30) 

(15.30) 
1 

(25.30) 
1 

(50.30) 

2.0 (0.6) 

6* 

- 0.2 

(1.6) 

2.7 

3.3 

- 0.3 

b 2 

(10,30) 

2.0 (0.6) 

6* 

- 0.2 

(1.6) 

2.7 

3.3 

- 0.3 

b 3 

(10,30) 

2.0 (0.6) 

6* 

- 0.2 

(1.6) 

2.7 

3.3 

- 0.3 







Table A1 

Coefficient Definitions for the Spectral Model (Al) - (A14) . 
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Table A2 

Definitions of the Fourier Coefficients Ai and Ti. 


n 

Ai i 

Ti 


\ f* U sin 2 (qz*)dz* f J* sin 2 (qz*)dz* 

- J r U* sin(qz* ) sin(nz* ) dz* \ \ J S- sin(qz*)sin(nz*)dz* 

X o ^ 0 
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3. RADIATIVELY-DRIVEN INTERACTIONS BETWEEN 
STRATOCUMULUS AND SYNOPTIC WAVES 


In a paper accepted for publication in the Journal of the 
Atmospheric Sciences (Clark, 1993), the effect of radiative 
cooling perturbations above stratocumulus on the structure of 
synoptic-scale waves is considered. The cloud position is linked 
to the waves via the phase of the low-level streamfunction. 
Coupling is strongest with cloud to the west of surface troughs. 
The resulting stationary structures capture the summertime 
pattern of stratocumulus off California and its linkage to the 
mid-Pacific ridge. The implication is that cloud-induced radiative 
cooling above cloudy regions plays an important role in 
determining observed structures over marine areas. 


A copy of this paper is attached. 
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ABSTRACT 


Quasi-geostrophic disturbances on a mid-latitude p-plane channel 
forced by radiative heating perturbations due to stratocumulus are 
considered. Longitudinal phase of the cloud is fixed to that of the low- 
level stream function. The background flow has a jet centered near the 
tropopause. Cloudiness is wavelike with cooling above cloudy areas 
that decreases exponentially with height. No perturbation cooling 
occurs above cloud-free areas. Forced steady waves are found. Though 
infinitesimal amplitude disturbances are considered, the problem is 
nonlinear because of the coupling between cloud and winds. The 
resulting structures are sensitive to the phase relationship between 
cloud and stream function with strongest coupling for cloud to the 
west of surface troughs. The waves have vertical scales roughly the 
troposphere depth. The stationary structures capture the summertime 
pattern of stratocumulus off California and its linkage to the mid- 
Pacific ridge. Zonal mean cloud cooling forces an adjustment to the 
mean westerlies that strengthens them in the northern half of the 
channel near the lower boundary. When this correction is allowed for, 
synoptic-scale amplitudes increase (decrease) just above the cloud in 
the northern (southern) part of the channel. Mean cloud cooling also 
renders the background potential vorticity gradient negative and thus 
baroclinically unstable just above the cloud in the northern domain. 
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1. Introduction 


Interactions between clouds and winds have important, but poorly 
understood, influences on atmospheric variability. Daily weather 
changes, as well as secular variations on inter-annual and climatic 
time scales, are modulated by such interactions. On the climatic scale, 
a major uncertainty in assessing the susceptibility of the general 
circulation to increasing CO 2 concentrations is the role of clouds. 
Lindzen (1990) doubts, because of the primitive state of our 
understanding of cloud-mediated feedbacks, whether it can be stated 
with certainty that global temperatures will increase in response to 
increasing CO 2 . Randall et al. (1984) show that cooling due to a 4% 
increase in global coverage by low-level stratiform cloud could more 
than offset the widely-predicted 2-4° K rise in global temperatures due 
to a doubling of CO 2 . 

There has been a lack of studies whose aim is to understand the 
interaction between motion systems and stratocumulus cloud. This is 
unfortunate because of widespread coverage by this cloud. For instance, 
according to Campana (1988), the largest contributor to total cloud 
cover at mid- and high latitudes is low cloud followed by middle and 
then high cloud. Also, Rossow and Schiffer (1991) show that in either 
hemisphere from 20° to 50° latitude most cloud tops occur at 
pressures of 800 mb or greater. 


Since the solar albedo of stratocumulus (typically 0.4) exceeds that 
of the underlying surface (typically < 0.1 for water and 0.2 for land), 
there is a deficit of absorbed solar radiation at the surface in cloudy 
regions. Near the ground, stratiform clouds also warm by absorbing and 
emitting infrared radiation. Ramanathan et al. (1989) show, however, 
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that the net radiative effect within the cloud and below is cooling. 

Thick decks of stratiform cloud also modify upwelling infrared 
radiation above the cloud. They emit less radiation than the lower 
boundary (especially over land areas) since their emission temperature 
is less than that of the underlying surface. The only exception is at high 
latitudes where cloud top temperatures can be higher than lower 
boundary temperatures and stratiform clouds can warm the 
atmospheric column above them. Randall et al. (1984) and Sohn and 
Smith (1992) show up to a 10% or 40 W m' 2 decrease in satellite- 
measured upwelling infrared radiances due to marine stratocumulus off 
the coast of California. 

Donner and Kuo (1984) calculated infrared heating rates allowing for 
low-level stratiform clouds and find, at mid-latitudes, a perturbation 
due to clouds of about 2 K d _1 near cloud top. They also find stationary 
quasi-geostrophic structures forced by topography, surface latent and 
sensible heating, and radiative heating (as modified by clouds). 
Resulting wave structures closely agree with tropospheric 
observations. Cloud induced asymmetries, they find, account for at 
least 20% of the calculated wave amplitudes. Donner and Kuo neglected 
meridional gradients of the mean zonal wind thus precluding wave 
energy accumulation in certain latitudinal belts because of meridional 
variations in refractive index (Karoly and Hoskins; 1982). This could 
lead to a significant error in their estimate of synoptic-scale 
sensitivity to cloud heating. There have been many other studies of the 
response of large-scale waves to stationary patterns of heating 
(Smagorinsky, 1953; Derome and Wiin Nielsen, 1971; Ashe, 1979; 
Hoskins and Karoly, 1981; Lindzen and Jacqumin, 1982; Aipert et al., 
1983). In these studies, clouds, if accounted for, are geographically 
fixed and not linked in any way to the disturbances they modify. 
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This study focuses on the effects of perturbations to the field of 
radiative heating due to stratiform clouds. The calculations are 
dynamical in that the cloud location is tied to the synoptic-scale 
disturbance it forces. Cooling rates due to stratocumulus will be very 
modest (a few deg. K d' 1 ). Calculated radiative cooling rates in very 
thin layers near the top of widespread cloud decks are roughly an order 
of magnitude larger (Nicholls, 1984). When this cooling spreads over a 
horizontal layer whose thickness matches the vertical scale of typical 

synoptic-scale disturbances, it is roughly of the magnitude used. 

* 

Certain stationary features of the large-scale cloud distribution 
link to synoptic scale wind patterns. For instance, persistent decks of 
marine stratocumulus occur in regions of subsidence to the east of 
subtropical highs at mid-latitudes (Schubert, 1976; Schubert et al., 
1979; Warren et al., 1988). In the northern hemisphere, this 
stratocumulus is most extensive during the summer when subsidence 
associated with the subtropical highs is strongest (Schubert et al., 
1979). In these regions, a strong temperature inversion with warm, dry 
air aloft caps the stratocumulus. Turbulence generated by shear in the 
cool moist air below the inversion and cloud top radiative cooling 
combine to maintain a sharp elevated base to the inversion layer (Lilly, 
1968; Schubert, 1976). Schubert et al. (1979) show that for a 
horizontally homogeneous steady-state cloud deck, the strength of the 
synoptic-scale divergence affects both cloud base and cloud top height: 
the stronger the divergence, the lower the cloud base and top and the 
thinner the cloud. Both sensible and latent heat flux divergences thus 
will be modified by divergence and thus there could be important 
feedbacks on the synoptic scale via this mechanism. Thus feedbacks of 
stratocumulus on synoptic scale motions are complex and are not solely 
due to low-level radiative cooling. This study concentrates on 
radiatively-driven coupling. 
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Other features of the stratiform cloud distribution closely couple to 
traveling weather systems. For example, the movement of cold air 
masses off the east coast of Asia and North America during winter 
behind propagating cyclones often leads to widespread marine 
stratocumulus. Also warm-frontal overrunning often creates broad 
regions of cloud ahead of propagating cyclones. 

It is thus reasonable to link the longitudinal phase of some features 
of the observed stratocumulus distribution to the phase of synoptic- 
scale wind and thermal patterns. The focus of this study is to examine 
the structure of synoptic-scale features that result from this coupling 
using linear quasi-geostrophic theory applied to a mid-latitude (3-plane. 
These structures will be compared with observations over marine 
regions where the effects of inhomogeneities in the underlying lower 
boundary are minimal. 

Assume the cloud distribution along a latitude circle is wavelike. 
Since cooling normally occurs in and above cloudy regions, mean and 
wavy components of cooling result, see Fig. 1. The mean cooling, if a 
function of latitude, can change the zonally-averaged winds and 
temperatures. They, in turn, can feed back on the synoptic-scale waves 
that couple with low-level clouds. This study will separately consider 
structures driven by mean and wavy components of the cloud cooling as 
well as feedbacks on stationary waves triggered by mean cloud cooling. 

Theory appears in Section 2 and results presented for an analytical 
constant background wind model. Stationary wave structures with a 
realistic background atmosphere are presented in Section 3. Forcing is 
by a combination of fixed lower boundary heating (to simulate 
geographically fixed patterns due to land-sea contrast) and cloud. 
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2. Theory 


The linearized quasi-geostrophic potential vorticity equation is 


^ + aQ 


aq 

at 


ax ay ax 



dq\ 


( 2 . 1 ) 


where \|/’ is the stream function. The vertical coordinate is log 
pressure, z=Hln(p 0 /p): H is scale height, p 0 = 1000 mb. The geostrophic 
winds Ug and Vg are, respectively, -3\|/73y and a\j f'/Bx. On a p-plane, the 
Coriolis parameter f=f 0 + py. The background density (p 0 ) is a function 
of z only. Linear dissipation in (2.1) arises from Newtonian cooling and 
Rayleigh friction. For convenience, the time scale of each is T<j = d'V 
The perturbation potential vorticity, q', relates to y' by 


q' 


v 2 y + 


Po az az ) 


(2.2) 


Define the heating rate, h' (dimension m s' 3 ), from the thermodynamic 
equation: 


+ n 2 w' = h' - db’ , 

9t ax ay dx 


(2.3) 
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where b' and B are, respectively, perturbation and background 
buoyancies. The buoyancy frequency squared is = dB/dz. The 
hydrostatic and thermal wind equations are, respectively, 

(2.4) 


* = » * . 

dz 


and 


au = _ J_3B 
dz f c By 


(2.5) 


The background potential vorticity gradient is 


3Q » 3fU + f|jL 

By dy 2 p 0 Bz 


Po au 
2 Bz 


\N 


(2.6) 


For steady heating at a single wavenumber k in the x direction, (2.1) 
and (2.2) combine to yield a diagnostic equation for y' ; 
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(2.7) 


where y’ and h’ are, respectively, 'P(y,z)exp(ikx) and n(y,z)exp(ikx) 
Dissipation is accounted for by letting 


U* = U(y,z) - i/(kT d ). 


( 2 . 8 ) 


Equation (2.7) is of the Helmholtz type and can usually be solved by 
successive over-relaxation. The local nature of its solution, neglecting 
dissipation for the moment, depends on the factor 
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In regions of the y-z plane where R is positive, a stationary solution is 
wavelike. This implies wave energy propagation. Where R is negative, 
energy propagation is not possible. 


a. Constant background wind 


Consider the response to heating where the wave and associated 
cloud have an arbitrary zonal phase speed c. Waves that neither grow 
nor decay with time, in spite of the radiative forcing, are considered. 
The wavy component of cloud heating is 


h' = n exp 


ik(x-ct) 


i^y - — 
D. 


The amplitude n is complex with modulus 2 K d -1 times g/0 o (g is 
gravity and 0 q = 273K). Heating decays exponentially with height with 
e-folding depth D=10 km but without any phase change in keeping with 
the model in Fig. 1. In this section and those to follow, only the region 
above cloud top, located at z=0, is resolved. The lower boundary 
condition of zero synoptic-scale vertical velocity is also applied at 
z=0. 

If the background wind is constant, dQ/dy = pin (2.6); (2.1) 
becomes 


where n 2 is assumed constant. If the perturbation stream function, y\ 
varies with height as exp(imz), the perturbation potential vorticity 
becomes q' = -(k 2 +^ 2 )\|/' - f 0 2 (m 2 +im/H)y7N 2 . 


In (2.10), y’ is the sum of a particular solution to the full equation 
(with heating) and a general solution without heating. If y’ has the 
horizontal structure exp(ikx+i^y), the homogeneous solution is B 
exp[i(kx+£y+mz)] provided 


m 2 


+ ini = N 
fo 


H a 



( 2 . 11 ) 


The constant B is not yet known. 

The particular solution is A exp[ik(x-ct)+i/?y-z/D] provided, using 

( 2 . 11 ), 


A = - 


n 



/ 

ikf 0 DU* 

m 2 + -3- + _3_ 

i h d 2 ch; 


( 2 . 12 ) 


The solution to (2.10) is now 

Y’ = exp[ik(x-ct)+tfy]|Aexpj--^J+ Bexp(imz)J. (2.13) 


The constant B follows from the thermodynamic equation applied at 
the lower boundary. Because of the log pressure vertical coordinate, the 



boundary condition of zero material vertical velocity must be applied 
with caution. The Appendix shows that a new non-Doppler term appears 
in the thermodynamic equation at z=0 for transient waves. It differs 
from the non-Doppler term discussed in Lindzen (1968) and Geisler and 
Dickinson (1975). They keep two terms arising from the expansion of 
material vertical velocity that cancel and thus obtain a different, but 
equivalent, formulation. The boundary condition at z=0 for a constant 
mean wind is 


fr 


i_ +u A|^.Mi?v + , nd ?v. h . 


\at dxj d z g 9 1 


3z 


(2.14) 


The non-Doppler term is the second on the left side. 


The constant B follows upon substituting (2.12) and (2.13) into 
(2.14): 


infm 2 +im+ 1 +-!&J 

1 H CH gDLTj 


J m 2 . im , 1 , 1 1 

imf 0 U* + 

f 0 N 2 c\ 

1 H D 2 CH/ 
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The stream function amplitude at the lower boundary is, using (2.12) 
and (2.15), 


A + B = - 


i n/m 2 + -Lin - Lm. + J-) 
\ h d m) 


\ 


kf ol m 2 + TL 


H 


D‘ 


2 + CHj 


= na (2.16) 


^-2- + imlT 

g 


The last equality defines G. Equation (2.16) enables the phase of the 
cloud field at z=0 relative to that of the stream function to be found as 
a function of U, k, d (or T<j), c, H, and D. Since the cloud is 180° out of 
phase with the heating (see Fig. 1), (2.16) states that phase{stream 



function} - phase{cloud} = arg{-G}. Suppose, for instance, arg{-G} - - 
90°. Cloud would be concentrated east of troughs and, if they tilt 
westward with height, on their warm side. Arg{-G} = +90° implies 
cloud to the west of surface troughs. Fig. 2a shows the sensitivity of 
arg{-G} to zonal wavenumber and background wind for stationary waves 
(c=0) with Td = 10 d, D=H=1 0 km and N=0.01 s' 1 . This value of T<j is in 
accord with the detailed calculations of Prinn (1977) for the radiative 
time scale of disturbances above an insulated lower boundary with a 
vertical wavelength about twice the troposphere depth. With westerly 
background flow, there is a broad region, especially for large 
wavenumbers, where the cloud mainly occurs to the east of the surface 
trough. Background easterlies or weak westerlies (especially for large 
east-west wavelengths or small zonal wavenumbers) favor cloud 
organized to the west of surface troughs. 

Wave stationarity has eliminated the non-Doppler term in (2.14). 
However, calculations show that even with traveling waves, this term 
is very small. Thus all solutions are mainly dependent on the Doppler- 
shifted background wind, U-c, and not independently on U and c. 
Therefore the ordinates in Figs. 2a and 2b could, without significant 
loss in accuracy, be U-c. 

The dimensionless surface stream function amplitude in Fig. 2b 
follows from (2.16) after dividing by g 2 /(f 0 N 2 ). Background Doppler- 
shifted westerlies and small zonal wavenumbers favor the largest 
response to stratocumulus cooling. One locus of maximum amplitude, 
labeled critical, occurs for U-c = 0. 

A second sloping locus of maximum response, labeled external 
Rossby, occurs in the region of westerly Doppler-shifted background 
winds. For zonal wavenumber five, according to Fig. 3, the vertical 


wavelength, 2n/Re(m), near this locus is over 100 km. The imaginary 
part of m, from Fig. 3, gives a 28 km e-folding depth for amplitude 
decay. Thus the mode is almost barotropic. If it was not for dissipation, 
the mean flow of about 8ms' 1 corresponding to this mode would 
render it as trapped or evanescent. Weak dissipation forces the wave to 
exhibit a slight westward tilt with height at all wavelengths and 
background winds. Stratocumulus primarily lies to the east of surface 
troughs for the external Rossby mode. Phase of the cloud field relative 
to that of the surface stream function is very sensitive to the Doppler- 
shifted background wind as witnessed by the concentration of constant 
phase lines near this locus in Fig. 4a. 

Marine stratocumulus off California lies to the east of the 
subtropical Pacific high in a region of weak westerlies. The pattern 
roughly corresponds to zonal wavenumber three to five and roughly 
conforms to the stationary critical mode structure found above. 
According to the theory, traveling stratocumulus patterns that 
radiatively couple to synoptic waves would take on the external Rossby 
mode structure with cloud concentrated to the east of surface troughs. 
Unfortunately there have not been enough observational studies of 
cyclonic systems to confirm or deny this linkage. 

b. Synoptic-scale vertical motion 

Summertime marine stratocumulus off the west coast of North 
America occurs in a region of synoptic-scale subsidence. The 
relationship between cloud location (as determined by the pattern of 
radiative cooling) and lower troposphere synoptic-scale vertical 
motion is implicit with the present analytical model. If the cloud 
coincides with descent, then just above the cloud radiational cooling is 
offset somewhat by compressional heating. Also, the capping inversion 



should be stronger in regions of synoptic-scale subsidence. 
Alternatively, if rising motion occurs in the region of low-level 
cloudiness, adiabatic cooling supplements cloud cooling. 


Substitute (2.14) into (2.3). The result is 


NV = 


IIe’ z/D - ikf 0 Ui- A- e’ z/D + irnBe imz 

D 


,ik(x-ct)+it y (2.17) 


After (2.12) and (2.15) are used for A and B, respectively, and (2.11) 
solved for m, w' can be found for a given heating rate amplitude n. 

Write the result of these substitutions into (2.17) as N 2 w' = Rh\ where 
R is complex and h'=nexp[ik(x-ct)+i^y]. The sum of cloud radiative and 

compressional heating is thus h'-N 2 w' = (1-R)h. Thus 1-R is the ratio 
of total heating to cloud radiative heating. 


The phase of 1-R at 1 km for stationary waves appears in Fig. 4 as a 
function of background wind and zonal wavenumber. For the reasons 
stated above, the background flow can, with little loss in accuracy, be 
considered the Doppler-shifted flow (U-c). In the shaded region, rising 
motion mainly occurs above cloudy regions that cool radiatively. 
Wavenumbers less than about six exhibit this anomalous behavior 
provided U-c < -5 m s . For U-c greater than about -5 m s and/or 
wavenumbers greater than about seven, low-level subsidence occurs 
predominantly in cloudy regions that cool radiatively. 


The latter pattern is most commonly observed in regions of 
persistent stratocumulus such as off the west coast of North America. 
Subsidence reinforces the capping inversion in these regions leading to 
a more permanent cloud deck. Again this simple theory agrees fairly 
well with observations of marine stratocumulus. 


1 4 



3. 


Influence of Background Wind Shear 


a. Response to wavy heating 

Consider steady waves forced by wavy heating on a mid-latitude 13- 
plane channel with a realistic background flow: 

U(y,z) = U 0 + U s cos(fy)exp[-(z-z 0 ) 2 /Du] , (3.1) 

with U 0 = 5 m s' 1 , U s = 15 m s' 1 , z 0 = 10 km, Dy = 10 km for z < z 0 and 

Djj = 20 km for z > z 0 , £=7t/3500 km' 1 . A 20 m s' 1 jet, Fig. 5a, centered 

at the tropopause and mid-channel, y=0, is described. The associated 
potential vorticity gradient, 3Q/3y, in Fig. 5b is positive everywhere 

with a maximum 2 km below the jet core. A stability jump occurs at 

the tropopause such that N=10'^ s' 1 (2.5x1 0'^ s' 1 ) for z < z 0 (z > z 0 ). 

The wavy heating has two components. 

• Fixed Pattern due to Land-Sea Temperature Contrast 

It is proportional to cos(kx)cos(^y) with amplitude 1 2 K d' 1 times 
g/0 o and occurs at the lower boundary. This heating is introduced by 
the thermodynamics Eq. (2.3) applied at z=0. Discretized Eq. (2.7) is 
solved in the interior on a 0.5 x 175 km height-latitude. The fixed 
heating does not contribute to the forcing term in (2.7) at any 
interior location. Thus this heating has a Dirac delta function 

1 An amplitude of 2 K d' 1 for the wavy component is equivalent to a cooling rate of nl 2 
times 2 K d' 1 or 3.14 K d _1 in the cloud/ areas and zero cooling in cloud free areas 
according to Fourier's theorem. 



distribution, 5(z), in the vertical. 


• Variable Phase Component due to Stratocumulus 

It links to the stream function at the lower boundary by an arbitrary 
phase O c ; if the stream function varies as cos(kx), the cloud cover 
varies as cos(kx-O c ). Negative cloud, of course, does not exist. Since 
cloud is 180° out of phase with radiative heating, see Fig. 1, O c = 

90° (-90°) signifies that stratocumulus coincides with low-level 
cold (warm) air advection to the west (east) of surface troughs 
assuming they tilt westward with increasing height. Thus <1> C equals 
arg{-G} discussed in Section 2a. Heating decays as exp(-z/D) where 
D = 10 km. Its formal representation is h'=IIexp(ikx) where 

n = n 0 , ¥ | Z= ^| expf-iOc - j-| cos(ly). (3.2) 

|\|/ (z=0)| l D) 

At mid-channel y = 0. Also n 0 = 2 K d'^ times g/0 o . 

Substitute (3.2) into (2.7); the steady response to steady wavy 
heating follows. Impose a radiation condition at the domain top (20 
km) and set the stream function equal to zero at the northern and 
southern walls, which are W=3500 km apart. 

There are some important points to note. 

• The fixed component of wavy heating is indispensable since it 
creates the disturbance to which the cloud and its cooling pattern is 
attached. That disturbance, in turn, is modulated in amplitude and 
phase by the cloud. 



• Although (2.7) is linear, condition (3.2) linking the cloud cooling and 
stream function renders the problem nonlinear. Sequential 
relaxation, nevertheless, was used to solve (2.7). Convergence, when 
achieved, was slow because cloud cooling was repeatedly adjusted 
according to (3.2) at each latitude. The algorithm succeeded for 
zonal wavenumbers greater than three. Planetary-scale solutions 
(wavenumbers one, two and three) were not found since these 
wavelengths are propagating (as will be shown later) at almost all 
channel locations. 

• In contrast to the constant background wind model, only the forced 
response to heating was found. It was not supplemented by a 
solution of the homogeneous problem. 

Wavy heating could trigger instabilities of the background state. 
Since, the potential vorticity gradient, dQ/dy, is positive everywhere, 
internal baroclinic instability (Charney and Stern, 1962) is precluded. 
Vertical shear near z = 0 is equivalent to a sheet close to the boundary 
where the potential vorticity increases to the south, Bretherton (1966). 
External baroclinic instabilities (James and Hoskins, 1985) are thus 
possible. They would have largest amplitude near the lower boundary 
and thus could have important implications on cloud-induced feedbacks. 
This study focuses on the steady response to heating and ignores 
spontaneous background instabilities. 

WKB theory is helpful in understanding solutions of (2.7). Karoly and 
Hoskins (1982) apply ray theory to interpret planetary wave 
propagation in regions of slowly (for WKB theory to be valid) varying R, 
defined in (2.9). Assume R is positive. Wave energy propagates along a 
ray or characteristic with a speed equal to the local group velocity. 
Rays tend to be parallel to the vector gradient of R. Wave energy thus 


focuses toward regions of maximum R. Zero wind lines where U = 0 and 
thus R is infinite will be locations of energy accumulation by 
stationary waves. As energy approaches the zero wind line, the group 
velocity approaches zero and the energy takes an infinite time to get 
there (at least according to linear, inviscid theory). 

Note, first, that locally unforced (n=0) solutions to (2.7) exist. In 
the inviscid WKB limit, the condition 


K 2 = i|?. JiL+jJ 

U 3V N 2\ 4H’j 


(3.3) 


must be satisfied, where k is the horizontal wavenumber, K^=k^-+£^, anc j 


p = m + -i— . (3.4) 

2H 

The stream function amplitude is wavelike in y and z: « exp(i^y+imz). 

The vertical wavenumber m can be complex, but p 2 is real. For p 2 >0, X P 
amplifies with height as exp(z/2H). An unforced solution is locally 

propagating (evanescent) if p^> 0 (p 2 <0). 

Equation (3.3) relates the horizontal and vertical wavenumbers of a 
stationary free wave. Let k cr be the zonal wavenumber at mid-channel 

(latitude 45°N) where the transition between evanescent (k<k cr ) and 
vertically propagating (k>k cr ) behavior occurs. Set p equal to zero in 
(3.3); then 


^2 1_ 3Q _ fp _ p2 

U 3y 4N 2 H 2 
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(3.5) 


Fig. 6 shows the integral zonal wavenumber n cr = rk cr cos(45°) for 
north-south wavelengths: Xy = W and 2W ( Xy=2K/C ). The earth’s radius is 
r. Mid-channel values of U and 3Q/3y are used. The jump in n cr across 
the tropopause reflects the sudden change in N. Ultra-long waves with 
wavenumbers one through four are propagating in most of the 
troposphere. The exception is very close to the ground where long-wave 
trapping occurs for Xy = W due to a maximum in curvature of U with 
respect to z. Short wave trapping for wavenumbers greater than six 
occurs for all z. All wavenumbers (except one) become or are trapped in 
the stratosphere. Turning points for longer wavelengths develop at 
larger altitudes than for shorter wavelengths. 

Consider a wave forced by stationary wavy heating, nexp(ikx), at 
z=0. Neglect dissipation. Use (2.5) for 3B/3y and apply (2.3) at z=0 
where w=0: 


ikfr 


I 


imU - 


V 


au 

dz. 


¥ = n. 


(3.6) 


Since the wave is stationary, w=0 is equivalent to setting the material 
vertical velocity equal to zero at the lower boundary. A resonant 
response to heating occurs if the bracketed term on the left side of 
(3.6) is zero. 


The mid-channel amplitude and phase of 1/(imU - 3U/3z) appear in 
Fig. 7 for Xy = 2W = 7000 km. The vertical wavenumber, m, follows 
from (2.11); acceptable roots have lm(|i)>0 thus preventing growth 
with z faster than exp(z/2H). The amplitude in Fig. 7 exhibits a broad 
peak near wavenumber five. In the absence of dissipation, the zonal 
wavenumber for a resonant response to the wavy heating is 5.35. The 
phase in Fig. 7 undergoes a shift of 100 to 110° near the peak response. 
Held (1983) shows that the resonant response is an external Rossby 


wave. With no dissipation, the resonant wave grows with z as exp(zU' 
"*3U/3z) near the lower boundary. In the mid- and upper troposphere, 
according to Fig. 6, the quasi-resonant response (wavenumber five) is 
propagating. A turning point occurs at the tropopause and the response 
is trapped in the stratosphere. These findings are consistent with 
those of Held (1983). 

Fig. 8 compares normalized amplitudes zonal wavenumber five for 
O c = -90° and +90°. Wave structures are in accord with deductions 
concerning locally evanescent and propagating behavior that follow 
from Fig. 6. Both structures grow in amplitude from z=0 to a maximum 
near the westerly jet core (mainly due to the decrease of background 
density). Wave evanescence forces decay in the stratospheric 
westerlies. The only significant difference between Fig. 8a and 8b 
occurs near domain base at mid-channel. With cloud to the west of 
surface troughs (O c = 90°), wavenumber five here is almost equal in 
amplitude to the jet core maximum. The O c = -90° structure is much 
more sharply peaked near the jet maximum; amplitudes near domain 
base at mid-channel are smaller than for G> c = 90°. The O c = 90° lower 
troposphere structure is amplified by cloud radiative cooling since cold 
air lying to the west of surface troughs (the waves tilt westward with 
increasing height) is cooled radiatively. The <X> C = -90° structure is 
weakened by radiative cooling since warm air located immediately 
above cloudy areas east of surface troughs is cooled. 

The Eliassen-Palm (EP) flux, shown in Fig. 9 for G> c = -90°, is a 
useful diagnostic of quasi-geostrophic disturbances (Edmon et al. 
,1980). Not only does the EP flux permit an evaluation of wave-mean 
flow interactions, it serves as a measure of wave propagation since it 
is parallel to the local group velocity vector. Provided WKB theory 
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holds, the group velocity is parallel to local wave energy flux vector. 
The EP flux is 


F = -UgVg j+ -k- Vgb k, (3.7) 

N 2 

/N /N 

where j and k are, respectively, unit vectors in the y and z directions. 
The overbar is a x-average over a wavelength. 

Streamline spacing in Fig. 9 is unrelated to flux magnitude; the 
arrows show flux direction and magnitude. Although the streamlines 
intersect the vertical walls 3500 km apart, there is no energy flux 
across them. Upward energy fluxes occur throughout the troposphere. 
Near mid-channel above 3 km, the flux vector is just about vertical. 
This implies, from (3.7), the dominance of the heat transport over the 
momentum transport contribution. The stationary wave must tilt 
westward with height. Because the regions of warm air advection 
ahead of surface troughs are cloudy, radiative cooling will decrease 
northward wave heat transport thereby making the EP vectors slightly 
less vertical than without cloud especially in the lower troposphere. 
Stratospheric EP fluxes are very small in Fig. 9 and mainly horizontal 
thus indicating wave evanescence. 

Sensitivity of stationary zonal wavenumber five to O c is examined 
in Fig. 10. Without cloud, forcing is due to land-sea contrast; the 
resulting mid-channel amplitude of v g at the lower boundary is, as 

shown, 14.7 m s'^. The phase without cloud follows from (2.3). 
Provided dissipation is small and wavenumber five is evanescent near 
z=0, stream function phase relative to that of the heating is, as 
displayed, -90°. Thus cooling (warming) due to land-sea temperature 
contrast lies to the west (east) of surface troughs. A balance is struck 
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between this heating pattern and advection by perturbation and 
background wind as follows: 

• Meridional advection of the background thermal field (background 
shear is westerly) gives further cooling (warming) to the west 
(east) of troughs. 

• Advection by the background westerlies gives warming (cooling) to 
the west (east) of troughs since troughs (ridges) are cold (warm). 

The surface wind with cloud, according to Fig. 10, varies by as much 
as 3 m s' 1 or 20% and surface trough or ridge locations by as much as 
14° longitude. The surface wave amplitude is larger with cloud to the 
west of surface troughs since cooling amplifies the perturbation 
thermal and geostrophic wind patterns. 


b. Response to mean cloud cooling 


There is a non-zero mean over a wavelength, n(y,z), of the rectified 
cloud cooling pattern in Fig. 1. Heating due to land-sea contrast, 
however, averages to zero. The response to mean cloud cooling follows 
from the zonally-averaged potential vorticity equation: 


d 2 W 


(po ds-L 

( i dQ 

dy 2 

Po dz 

In 2 dz l 

lu* d y 



if o 3 Po n | 

kll’po 3z l N 2 ) 


(3.8) 


The stream function, V F, gives the mean zonal wind and buoyancy 
according to - ST 1 /dy and foB'P/^z respectively. Dissipation of potential 
vorticity in (3.8) balances cloud radiative generation. 


The boundary condition at z=0 follows from the zonally-averaged 
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thermodynamic equation: -df 0 9'*'/3z + n = 0. Mean cloud cooling is small 
enough at the upper boundary that the zonal mean buoyancy perturbation 
or d'V/dz is zero. At lateral domain boundaries, the zonally-averaged 
ageostrophic meridional flow must be zero. Thus the only effect left in 
the zonally-averaged steady u momentum equation is Rayleigh friction. 
The lateral boundary condition is thus 34 / /3y = 0. 

Sequential relaxation solves (3.8). Resulting tropospheric zonal 
wind, buoyancy, and vertical velocity fields appear in Fig. 11; the latter 
follows from the zonally-averaged thermodynamic equation. Mean cloud 
cooling creates a negative buoyancy perturbation, Fig. 11b, that is in 
thermal wind balance with a zonal wind whose maximum speed is 2.7 
m s’ 1 at the lower boundary, Fig. 11a. The wind is westerly in the 
northern part of the domain and easterly to the south. In most of the 
channel, compressional heating in weak subsidence of up to 0.06 cm s' 
1, Fig. 11c, offsets cloud cooling. Rising motion of up to 0.09 cm s' ^ 
occurs near the northern and southern boundaries. Since there is no 
cloud cooling at these boundaries, adiabatic cooling is balanced by 
Newtonian warming. 

c. Feedbacks of mean cloud cooling 

Linearization decouples the responses to wavy and mean cloud 
cooling since they are assumed to be of the same order of a 
dimensionless amplitude. Thus the stationary wave forced by wavy 
cloud cooling does not sense the mean flow, and accompanying thermal 
field, set up mean cloud cooling. Useful insights are gained by departing 
from strict linearity to find the stationary wave response to 
background flow alterations by mean cloud cooling. 

Add the zonal flow forced by mean cooling, Fig. 11a, to the 
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background flow of Fig. 5. The jet core, Fig. 12a, shifts southward and 
westerly vertical shear just above cloud decreases (increases) in the 
northern (southern) half of the channel. The result, near the northern 
boundary, is little vertical shear for all z. Substantial changes in 
potential vorticity gradient, 9Q/9y, occur near z=0 (see Fig. 12b). A 
region of negative gradient materializes in the northern half of the 
domain. Thus mean cloud cooling has rendered the background flow 
baroclinically unstable (Charney and Stern, 1962). At the same time, 
BQ/dy doubles immediately above the cloud to the south. 

The difference between stationary wavenumber five amplitude 
obtained with modified and unmodified mean flow appears in Fig. 13. 
Mean cloud cooling affects the structure below 5 km. Wave 
amplification (decay) occurs in the northern (southern) half of the 
channel. The region negative dQ/dy to the north strongly influences the 
local wave amplitude. According to (3.5), all wavenumbers are trapped 
when 3Q/3y<0 provided U>0. Thus wave amplification occurs above 
cloud top as the wave cannot propagate energy laterally or vertically 
away from the energy source as efficiently as with the unmodified 
flow. To the south, according to (3.5), the critical wavenumber, k cr , for 
the transition from wave propagation to evanescence increases in 
response to doubling 3Q/9y. Propagation of energy away from the energy 
source due to wavy cloud cooling is facilitated. The local wave 
amplitude thus decreases. 
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4. Discussion 


Significant progress has occurred recently in understanding the 
mid-latitude planetary boundary layer (PBL) when cloud is present. 
Most studies of the cloud-topped PBL treat synoptic scale wind, 
divergence, and advection as given; they focus on the ensuing PBL 
thermodynamics of vertical mixing under horizontally homogeneous 
conditions (Lilly, 1968; Schubert, 1976 and 1979; Nicholls, 1984). 
Because of PBL vertical mixing (of heat, momentum and water 
substance), mass entrainment across the capping inversion, and latent 
heating through drizzle formation, the PBL can in turn influence the 
synoptic mass, pressure and wind fields. The present model crudely 
attempts to close this feedback loop. Without explicit moisture and a 
PBL, this study focuses on the feedbacks to the synoptic scale of a 
prescribed pattern of low-level stratiform cloud that is linked by 
longitudinal phase with the in situ synoptic wind field. The 
hypothesized process responsible for this linkage is cloud-induced 
perturbations to the field of radiative heating. 

It is reasonable to not account for the PBL structure if the vertical 
scale of the disturbances excited by cloud cooling is much larger than 
the boundary layer thickness. For the linear structures forced by wavy 
cloud cooling whose vertical scales are on the order of the troposphere 
depth, this is a posteriori a reasonable assumption. Infrared cooling 
due to stratocumulus, however, is sharply peaked near cloud top. A 
caveat that can be attached to the present study is that a considerably 
smaller e-folding depth for the decay of cloud cooling than 10 km 
should have been used. Wavy cloud cooling with a 2 km e-folding depth 
forces stationary synoptic waves whose vertical and meridional 
structures are almost identical with those in Fig. 8. The conclusion is 
that the background zonal wind field plays a dominant role in fixing 
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linear wave structures. The response to mean cloud cooling with a 2 
km e-folding depth, however, is considerably shallower and it is hard 
to justify ignoring the PBL structure. It is not within the scope of this 
study to pursue this matter any further. Suffice it to say that the 
rectified nature of cloud cooling associated with wavy stratocumulus 
patterns could have implications that go far beyond those pointed out in 
this study. 

A second caveat is that the cloud top radiative cooling perturbation 
used in the present calculations was very modest (a few deg. K d""'). As 
pointed out in the introduction, cooling near cloud top can be an order 
of magnitude or more larger (Nicholls, 1984). It is thus possible that 
the response to cloud cooling is highly nonlinear near these 
concentrated layers of cooling. The resulting synoptic waves would be 
strongly coupled to the PBL and explicit representation of the PBL 
should be necessary to fully account for their structure. 

The main finding of this study is that the observed coupling between 
summertime stratocumulus and northern Pacific surface winds and 
pressures is plausibly driven by cloud radiative cooling perturbations. 
The model produces stationary structures that resemble observed 
patterns not only in the phase between cloud and low-level horizontal 
winds but in the relationship between low-level subsidence and cloud. 

There are other mechanisms besides radiative cooling that 
influence the linkage of the cloud-topped marine PBL with synoptic 
disturbances. For instance, Schubert (1976) showed that synoptic 
divergence influences PBL depth in a steady mixed layer PBL model. 
Horizontal divergence, presumably, would increase turbulent sensible 
and latent heat flux divergences by lowering the PBL top. A feedback on 
the synoptic scale motion could ensue. Another mechanism could be 
drizzle formation. Linder conditions of strong cold air advection, drizzle 
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often forms in stratocumulus (Nicholls, 1987). A net latent heating of 
the PBL follows while, at the same time, the PBL is stabilized through 
heating in the cloud layer and evaporative cooling below. Sea surface 
temperatures also undoubtedly influence observed stratocumulus 
patterns. 

It is unlikely, however, that these factors, collectively or 
individually, can account for the observed coupling over the 
summertime northern Pacific. Cloud radiative cooling must be 
accounted for in any comprehensive theory that explains the observed 
linkages. 
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APPENDIX 


The material vertical velocity (w*) must be zero on a horizontal 
solid boundary. The vertical velocity (w) in the log pressure coordinate 
system is 


w = - — co , (A1 ) 

P 

where to = dp/dt. If geometric height, z*, is the vertical coordinate, co 
can be expanded as 


0> = i*£ + Vg-Vp + W'^2-, (A2) 

at 9 3 Z * 

where, since the horizontal wind is geostrophic, the middle term on the 
right side drops out. At the lower boundary w* is zero. Thus, from (A1) 
and (A2), 


w--tL^. (A3) 

p at 


In linear theory (A3) applies at z=0 as well as at z*=0. The quasi- 
geostrophic stream function \|/ relates to pressure by \jr=p/(f 0 Po) » ^ us 
(A3) becomes 


w = .H = J±?V 

v at Vo at ' 


(A4) 
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The last equality follows after linearization and defining vy 0 as a 
background streamfunction (=Po/po^o)- After the scale height (H) is 
introduced, then at the lower boundary, where w* is zero, 


w = 


9 at ' 


(A5) 


Thus, for example, the zonally-averaged thermodynamic equation at the 
lower boundary is 


aB f 0 N 2 ay _ h 

at g at 


(A6) 
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FIGURE CAPTIONS 


Fig. 1 Schematic of wave structure showing location of cloud 
cooling. 

Fig. 2 Phase (degrees), (a), and amplitude (dimensionless) in units of 
10 6 , (b), as functions of background wind (U) and zonal 
wavenumber of stream function at z=0 in analytical model. 
Dissipation time is 10 days; phase speed is zero. N=0.01 s' 1 . 
Shaded region in (b) is where amplitude is less than 16. 

Fig. 3 Real and imaginary parts of vertical wavenumber in units of 
10'® m' 1 for stationary waves of zonal wavenumber five in 
analytical model. N=0.01 s' 1 . 

Fig. 4 Phase of factor 1-R (defined in text) at z=1 km 

for zonal wavenumber five, = 10 d and U=10 m s' 1 . 
Meridional wavelength is 7000 km or 2W, where W is channel 
width. Phase is labelled in degrees. Shaded region is where 

phase is less than -90°. 

Fig. 5 (a) Background wind, U(y,z); interval 2.5 m s' 1 , (b) Potential 

vorticity gradient, 3Q/3y; interval 10' 11 m' 1 s' 1 . 

Fig. 6 Critical zonal wavenumber separating vertically propagating 
from evanescent solutions at mid-channel. Labels Xy=\N (2W) 
refer to wavelength in y direction where W is channel width. 

Fig. 7 Amplitude (s' 1 ), and phase, of 1/(imU-3U/3z) at lower 



boundary in mid-channel versus wavenumber. 


Fig. 8 Amplitude of perturbation stream function (normalized as 
maximum value is unity) for zonal wavenumber is five. 

(a) <£> c = -90°; (b) O c = 90°. 

i 

Fig. 9 Eliassen-Palm flux for wavenumber five and O c = -90°. 

Fig. 10 Perturbation geostrophic wind at lower boundary at mid- 
channel, Vg, and phase of stream function, X, at lower 
boundary as functions of the phase of stratocumulus-induced 
heating field. If heating due to land-sea contrast goes as 
cos(kx), stream function is cos(kx+>.). Also shown are Vg and X 
if clouds do not occur and the wave is driven by fixed surface 
heating only. 

Fig. 11 Tropospheric response to mean cloud cooling (a) Zonal flow 

forced by mean stratocumulus cooling. Contour interval is 0.6 
m s' 1 , (b) Buoyancy perturbation; interval 0.03 m s"^. (c) 

Vertical velocity; interval 0.03 cm s' 1 . Dashed (solid) 
contours indicate negative (positive) values. Zero contours are 
labeled. 

Fig. 12 (a) Background wind supplemented by contribution due to mean 

cloud cooling; interval 2.5 m s' 1 , (b) Potential vorticity 
gradient (negative region shaded) corresponding to flow in (a); 
interval 10' 11 m' 1 s' 1 . 

Fig. 13 Change in stream function amplitude when mean flow forced by 
cloud cooling is added to background flow; wavenumber five. 
d> c = 90°. Values are normalized to unity. 
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